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ABSTRACT 

Two dereverberation techniques are applied to synthetic 
seismic d^ta and their performance in removing water column 
multiples is evaluated and compared. One method is an applica- 
tion of homomorphic deconvolution and the other utilizes 
linear estimation based on a minimum mean square error criterion. 
The analytical formulations of both methods are discussed. 
Performance is evaluated in terms of three criteria; percent 
of multiple energy removed, percent of signal (reflector) 
distox'tion, and visual improvement of the data. Results are 
presented which represent the performance of both algorithms 
for a range of environmental and signal processing parameters 
including white noise level, multiple coherence, reflector/ 
multiple overlap, filter parameters and water column travel 
time estimate. The techniques are found to have comparable 
effectiveness on the synthetic data; however, indications are 
that homomorphic dereverberation has greater potential in 
shallow water applications while the linear technique appears 
to be more efficient for deep water data. 
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CHAPTER I 



INTRODUCTION AND PROBLEM STATEMENT 



The geologic structure of the earth beneath the seafloor 
is most often determined by seismic profiling. The procedure 
generally involves excitation of an impulsive acoustic source 
near the sea surface and recording of the reflected earth 
response with a hydrophone array. Low frequency sound 
penetrates the bottom and propagates in the substrata with 
reflections occurring at discontinuities in the acoustic 
impedance of the earth. The thickness and density of sub- 
bottom layers may be estimated from the reflected acoustic 
signal, or seismogram. 

The earth is modelled as a discrete layered medium with 
distinct interfaces for most seismic applications. This 
assumed structure, while not strictly accurate, has led to 
good processing results in practical seismic work and has the 
additional advantage of being analytically tractable. We shall 
employ this assumption throughout the present analysis. A 
detailed description of the earth model used is given in 
Chapter III. 

The amount of energy reflected at a discontinuity is 
ideally measured by the reflection coefficient. 
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where and c^, C 2 are the densities and acoustic 

velocities in media 1 and 2, Here the sound propagates from 
medium 1 to medium 2 and normal incidence has been assumed. 

The assumption of normally incident plane waves is generally 
made in single channel (one receiver.) seismology since it 
leads to analytical simplicity and appears to be reasonably 
accurate. In utilizing this simplified model we. have ignored 
near field effects and spreading losses. These are not of 
major importance in single channel dereverberation and their 
inclusion would unnecessarily complicate the earth model. 

It is evident from the reflection coefficient expression 
that large reflections occur at points of significant change 
in impedance. The largest impedance discontinuity encountered 
by seismic signals occurs at the water-air interface where R 
is nearly -1. The water-bottom interface is also a strong 
reflector in most cases. Thus, the water column becomes a 
reverberating channel wherein a significant portion of the 
source energy is trapped. Repeated reflections from the bottom, 
or multiples, are received at intervals corresponding to the 
two-way travel time of sound in the water column. Deeper 
reflections from the substrata are masked by multiples when 
their arrivals are nearly coincident in time. Since propaga- 
tion losses are much smaller in water, the energy in the 
multiples is usually large compared to that in the deeper 
reflectors. Thus, water column multiples form an unwanted 
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component of the seismogram. 

The reflection coefficient expression indicates that 
all earth layers also introduce reverberation. Except in some 
shallow water situations, internal multiples are generally 
not a serious problem for two major reasons. Most importantly 
the acoustic attenuation in the earth is much greater than 
that in water so that very little internal multiple energy 
is actually returned to the receiver. Secondly, the reflection 
coefficients at earth layer boundaries are usually small 
compared to those at the surface and seafloor so that a 
relatively small part of the incident energy is actually 
trapped. 

Figure 1 shows a synthetic seismogram with strong 
multiples. Response amplitude is measured on the ordinate 
and travel time in seconds on the abscissa. The first large 
signal component is the bottom reflection at one second of 
travel time or about 750 meters v/ater depth. The first 
multiple is an attenuated, phase-inverted replica of this 
reflection at two seconds. Note that the return from a 
reflection horizon at two seconds travel time would coincide 
with this multiple and be obscured. The overall periodicity 
of the multiples is apparent in this plot. Actual reflectors 
occur at 1.5, 2.2 and 2.9 seconds. Figure 2 shows the seismic 
environment which would produce such a seismogram. 

The first practical multiple analysis and dereverberation 
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Figure 1 Synthetic seismogram including three multiples. Reflectors are designated 
Ri and multiples by Bi. 
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Figure 2 Earth structure v;hich v>70uld produce seismogram 
like that shown in Figure 1. 
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algorithm was proposed by Backus [1]. His approach was to 
characterize the water column as a sharp, ringing filter v/ith 
an impulse response composed of a v^eighted sum of delayed 
impulses. The weights and delays are determined by the bottom 
reflection coefficient and the water column travel time 
respectively. This model leads to a three-point operator v;ith 
elements spaced at intervals corresponding to the two-way 
water travel time. Implementation of the Backus filter requires 
estimation of the bottom reflection coefficient and water 
column travel time. Several aspects of the performance of 
this method are discussed in Chapter II. 

Spatial processing has also been used to reduce 
multiples [2]. Spatial schemes normally require multichannel 
arrays of large physical extent which can be effectively focused 
to discriminate against reverberation. Such systems are widely 
used and quite effective in shallow water but their costs, both 
for hardv;are and data processing, are very high. Hence, there 
is still a need for time domain multiple removal techniques in 
deep water situations and for single channel systems. 

Two techniques have recently been applied to seismic 
multiple removal with demonstrated success. The first, an 
inverse filter algorithm based on a tapped delay line model, is 
due to Baggeroer [3j, and is referred to hereafter as the TDL 
filter. A tapped delay line is simply a realization of the 
time domain convolution of a signal and a gapped operator [4]. 
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Results of this method have proven superior to those of three- 
point filtering, at least in some deep water situations. 

A nonlinear filtering scheme using homomorphic deconvolu- 
tion has been applied to several aspects of seismic signal 
processing. The use of this method for dereverberation has 
been dem.onstrated by Stoffa, Buhl and Bryan [5]. The homomorphic 
transformation is essentially a mapping from convolution to 
addition so that, after transforming, deconvolution can be 
accomplished by simple linear filtering. Seismic dereverbera- 
tion appears theoretically to be a very promising application 
of homomorphic deconvolution because of the distinct properties 
of seismic signal components. The method has not, however, 
been fully evaluated or widely used in practice. 

The motivation for this study arises from the disparate 
theoretical mechanisms by v;hich these techniques operate to 
perform the same function. Since analytical comparison is 
not feasible, this functional, comparative approach is 
thought to be the best means of gaining insight into this 
interesting problem. 

The purpose of the analysis is twofold. First, it is 
intended to indicate those factors which have significant 
effects on the performance of each algorithm. The factors 
to be considered are environmental variables and processing 
parameters. These are discussed in Chapter III. Secondly, 
the analysis is intended to point out the relative strengths 
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and weaknesses of the methods by comparison of results on 
similar data. 

Each algorithm is evaluated for a range of simulated 
processing conditions. Quantitative and qualitative criteria 
are specified which provide a comprehensive description of 
the manner in which each signal is affected by processing 
for multiple removal. These criteria also serve, as a basis 
for comparison of results. The scope of the analysis and 
the specific performance criteria are discussed thoroughly 
in Chapter III. 
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CHAPTER II 

ANALYTICAL FORMULATION AND IMPLEMENTATION 

A . Analytical Formulation of the TDL Filter 

A simple feedback model for the propagation of reverberat- 
ing seismic signals is given by Baggeroer [3] . Multiple 
removal based on this model is then formulated as an inverse 
filtering problem. The dereverberation filter is designed 
using a least squared error criterion and the constraint that 
the filter have a tapped delay line structure. 

The formulation will be developed here from a different 
point of view using Baggeroer 's feedback model as a starting 
point. A summary of the feedback model is included for clarity. 

Figure 3 shows the Laplace transform representation of a 
propagating seismic signal. S(s) is the transform of the 
source signature. The signal first encounters the downward 
travel time delay which corresponds to a phase shift in this 
domain, represents the transfer function of the earth 

beneath the water column including the reflections from layer 
boundaries which comprise the desired information. Internal 
multiples or reverberations between the various earth layers 
are also included in Hj^(s). Another phase shift corresponds 
the return of the reflected signal through the water column. 

P(s) is the feedback gain representing the water column multiple 
mechanism. In most cases this is well approximated by -1, 
which corresponds to the nearly perfect pressure release 
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Figure 3 Laplace transform model of propagating seismic signal 
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reflection at the surface. ^^qCs) includes the observation 
effects such as hydrophone bandwidth and array tov/ depth. 
Ambient noise and reciever front end noise are modelled as 
additive white Gaussian noise. 

The overall transfer function is 

-2sT 



R^(s) H^(s) U^is) e 

S(s) 1 - P(s) Hj^(s) e“^^w 



( 1 ) 



where T is the one-way water travel time and R (s) is the 
w o 

received signal without additive noise. 

It is apparent that the presence of multiples is due only 
to the denominator of this expression. Thus far v;e have 
assumed implicitly that the earth response can be modelled 
accurately as a linear system and that the multiples are 
exactly periodic. The validity of these assumptions will 
become apparent in the discussion of performance in Chapter IV. 

The obvious task is now to design an inverse filter having 
the form 

“2sT^ 

F(s) = tl - P(s) Hj^(s) e . 

Hence, we are required to estimate T^ and the impulse response, 
hj^(t), corresponding to Hj^(s). The earth response need not 
be estimated precisely for its entire duration. Estimating 
the dominant energy part of hj^(t) is adequate to produce an 
effective dereverberation filter. A typical deep water 
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seismogram has the great majority of its energy concentrated 
in the first 200-300 msec of its duration. Effective represent- 
ation of this portion of the signal requires about 10-20 
filter coefficients, depending on the bottom and source 
characteristics . 

The transfer function of equation (1) can be re-written 
in series form as 



R (s) “ ,, -2nsT 

-2 = H (s) I (-1)^ ^ H (s) e . 

S(s) ° n=l 

The received signal then has the time domain representation 



r (t) = s(t) * [h^(t) * I hj^(t - 2nT^,)] 

n=l 



where * represents convolution. This can be rewritten as the 
sum of the primary return and the multiple signal. 



ro(t) = 



s (t) 



* ho(t) 









n=2 



'] 



r^(t) = b(t) + m(t) 
where 



b(t) = s(t) * h (t) * h. (t-2T ) 

o b w 

is the received primary and 

00 

m(t) = s(t) * h (t) * y h, (t-2nT ) 

o ^„b w 

n=2 
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is the received multiple r.ignal. 

The estimation problem, given the feedback model of 
figure 3, is to determine b(t) in the presence of m. (t) and 
white noise, v;Ct). It is convenient to group the unwanted 
signal components. 



Since the unwanted component is an additive one, we 
can consider estimating n(t) and subtracting the result from 
the received signal. We then have the filtering problem 
depicted in figure 4. 



Here f(t) is the filter impulse response and n(t) is the 
minimum m.ean square error (Mt^SE) estimate of n(t), given the 
received signal r(t). 

The number of digital filter coefficients to be estimated 
is 2TI'7+1, where T is the effective duration (portion contain- 
ing about 80% or more of the signal energy) of hj^(t) and W 
is the signal bandwidth. The coefficients v/ill then be 
spaced at the Myquist sampling interval of 1/2W seconds. 



n ( t ) = m ( t ) + V7 ( t ) 



( 2 ) 



r(t) = b(t) + n{ 




Figure 4 
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The optimum digital filter for n(t) v/ill have coefficients, 
f , which minimize 



e = E { I 
i=j 



n(i/2W) - n(i/2W) 



(3a) 



v/nere 



n(i/2W) = I f,.*r[-^-^^^ - 2 t 



k=i 



2W 



v; 



(3b) 



The input in (3b) is shifted by the two-way travel time to avoid 
useless filtering of the signal prior to the bottom reflection. 
[i^/2W, i^/2V7] 3S the time interval over v/hich n(t) is observed. 



Substituting (3b) into (3a) yields 



e = E 



f r 

I I I 

i i=i L k=i 
o o 






(i-k) 

2W 



- 2T 



w 



- n(i/2W) 



Minimizing , 
9e 



= E 



8f 



r f 

[m. 

1=1 






(i-k) 



2V1 



-n (i/2W) 



-r 



“ (k/2W+2T„) . 



(i-k) 

2W 



(4) 



“ 2T. 



w 



Here we have assumed stationarity over the duration of the 
multiple period. This assumption has led to effective process- 
ing of both real and synthetic data. From (2) 



R (k/2W+2T ) = R (k/2W+2T ) + R (k/2W+2T ) 
nr ^ w mr ' w wr w 
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Here it is useful to assxune (see ref.£3j, p.l5) that for shifts 
close to 2T^ the cross-correlation of w(t) and r(t) is very 
small compared to which will generally have a peak in 

this region. This is equivalent to assuming that the v^hite 
noise has a very short correlation time compared to 2T^^. We then 
have 



R (k/2V7+2T R (k/2W+2T ) 

nr w mr w 

so that (4) becomes 
^f 

I fy. R^^((i-k)/2w) = R^^(k/2W+2T ) (4a) 

k=i 

o 

Baggeroer has derived equation (4a) by designing the 
Wiener filter for b(t) with the constraint that the filter have 
a tapped delay line structure, i.e. 

2TW 

f(t) = 6(t) - ^ fy. 6 (t-k/2W-2T^) (5) 

k=o 

I'Then our estimate, n(t) , is subtracted from the unshifted, 
received signal, the resulting overall filter operation has 
exactly the form of (5) • This indirect approach yields the 
estimator equations of reference [3]’ without imposing the TDL 
structure directly. 

The above derivation also emphasizes the estimator- 
subtractor or prediction error structure of this filter. The 
entire impulse response may be written as follows : 
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2T zeros 
w 



This is the prediction error structure for a prediction distance 

of 2T^. Equation (4a), however, which generates the 

differs from the prediction equations in that the right hand 

side vector is R (t+ 2T ) rather than R (x+2T ). As written, 

mr w rr w ' 

equation (4a) corresponds to the Wiener filter which produces 
the MT4SE estimate of m(t) with r(t-2T ) as an input. Subtraction 
of this estimate from r(t) whitens only those spectral compon- 
ents. which are due to the multiple. 

It is simply proved that the magnitude of the error in b(t) 
is equal to that in the prediction operator. 

b(t) = r (t) - n (t) 

b(t) = b{t) + (n(t)-Mt)) 

|b(t)-b(t)l = |n(t)-n(t)[ 

Thus far the only departures from optimum estimation have 

been the two assumptions of stationarity • and the relative 

insignificance of R (t+ 2T ) . One further assumption is 

wr w 

required for actual implementation of the filter. Note that 
R (t+ 2T ) is a required input which is apparently not measur- 
able from the given data. Baggeroer has observed that for the 
deep water case, which is of primary interest for this method. 



I 
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B^r(T+2T^)^R„(T+2T^) . 

That is, for shifts of nearly twice the v/ater travel time, the 
great majority of the crosscorrelation energy is due to m(t) . 
Hence, the equations used for data processing are given by 

^f 

I fv R ((J-k)/2W) = R (j/2W+2T ) j=0,l,...2TW 

Having seen the analytical formulation of this inverse 
filtering procedure it is instructive to compare it with the 
Backus three-point method. Processing actual data with both 
filters (ref. [3]) has shown the Backus filter to be significantly 
less effective. Some reasons for this are apparent from the 
foregoing analysis . 

The Backus filter is rigidly dependent on the accuracy of 
two assumptions. The first, the assumption of strictly periodic 
multiples, is violated due to the horizontal separation of 
source and receiver. This effect becomes more severe as v/ater 
depth decreases. Since the Backus filter is implemented as 
only three, equidistant operator coefficients it is very sensitive 
to this lack of periodicity. Even if the statistics of the 
signal generate very accurate estimates of the bottom reflec- 
tion coefficient the filter structure is so simple and rigid 
that proper cancellation v/ill not occur if the multiples are 
significantly aperiodic. 
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A second restrictive assumption of the Backus filter is 
that the bottom reflection should be accurately characterized 
by a single reflection coefficient. All of the statistical 
information available is forced into a single parameter 
estimation scheme. It is apparent that such a filter lacks 
flexibility for dealing with more complicated bottom interaction 
mechanisms . 

The relatively better performance of the TDL filter on 
real data is apparently due to its greater inherent flexibility. 
That is, the finite length impulse response, or prediction 
operator, gives the filter a capability for removing reverbera- 
tion effectively in cases where the bottom response is not 
accurately modelled as a weighted impulse. If the bottom has 
a ringing or smearing effect on the incident signal then the 
deconvolution operator must be extended in time. The Backus 
filter, because of its rigid structure, cannot accomodate these 
situations. The TDL structure provides 2TW+1 (usually 10-20) 
parameters v/hich can be varied in the design procedure to 
optimize de reverberation of each seismogram. The special case 
of an ideal bottom will generate a filter response which is 
essentially a single spike proportional to the bottom reflection 
coefficient. This result has been confirmed in the analysis of 
synthetic data. In such a case the TDL filter consists 
basically of the first tv;o points of the Backus three-point 
filter. Performance (multiple energy removed) in these cases 
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was found to be essentially independent of operator length. 

The structural flexibility of the TDL filter also gives 
it some potential for dealing with aperiodic multiples. It 
should be noted, however, that the TDL algorithm, like the 
Backus and other classical dereverberation techniques, is 
essentially a correlation-cancellation operation. Consequently, 
increased aperiodicity of multiples can be expected to degrade 
performance in all cases. 

B. Implementation of the TDL Algorithm 

Figure 5 shows a flow diagram of the filter implementation 
used for this analysis. Actual programming was done in 
Fortran IV for use on a 32K computer. Referring to figure 5, 
the correlation function is computed by the standard shift-and- 
add operation with no windowing applied. Results of windov.’ing 
are included in reference [3]. Correlation time is variable 
and may be specified by the operator. The crosscorrelation 
function is approximated by the correlation function as discussed 
in the previous section. 

Solution of the filter equations is accomplished by 
conventional matrix inversion. The Toeplitz symmetry can be 
exploited for computational savings. Spacing of the operator 
elements is determined by the estimated signal bandwidth which 
is specified as an input parameter. 

Actual deconvolution is implemented exactly as shown in 
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Figure 5 Flov/ chart of TDL algorithm. 
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Figure 6 Tapped delay line model of convolution with a gapped operator. 
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figure 6, i.e., by means of a tapped delay line or, equivalently, 
convolution v;ith the gapped operator. 

C. Analytical Formulation of Homomorphic Deconvolution 

A homomorphic system is one which obeys a generalized 
principle of superposition and which can be represented as an 
algebraically linear transformation between two vector spaces. 

A detailed description of the theory of homomorphic systems 
is given by Oppenheim and Schafer [6]. This material will not 
be repeated in depth here; rather, we shall discuss the basic 
characteristics of homomorphic systems for convolution with 
emphasis on those properties which facilitate dereverberation. 
Additional discussions of these properties are found in 
references [5] , [7] and [8] . 

The usefulness of linear systems for' separating additively 
combined signals is due primarily to superposition. Signals 
which are added and happen to be disjoint in the frequency 
domain can be separated by means of an appropriate bandpass 
filter. Homom.orphic systems for convolution have a similar 
effect on signals which have been convolved. That is, a 
homomorphic transformation maps the input signal to a domain 
in which the convolved components may be disjoint. Such a 
transformation is illustrated in figure 7. 




Figure 7 
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The mapping characteristic of this system is intuitively 
apparent. Suppose x(n) is composed of two components v/hich 
have been convolved. 

X (n) = x^ (n) * X2 (n) 

The z-transform operation yields a signal with multiplied 
components, X^(z) and X2(z). The logarithm output is 

X(z) = log[X^(z)] + log [X2(z)]. 

The inverse z-transform, 

/N /N /\ 

X (n) = x^ (n) + X2 (n) , 

preserves the additive combination of the components and yields 

a sequence which is real and stable for a real and stable x(n) . 
/% 

The sequence x(n) is called the complex cepstrum of x(n). 
Although it is real for real inputs, the "complex" is retained 
to emphasize that it contains both the magnitude and phase 
information from X(z). Hence the complex logarithm is required 
even for real x(n). (We shall omit the modifier here for 
brevity. ) 

The cepstrum variable, n, is normally called the quefrency 
(a paraphrase of frequency) or period variable. Filter opera- 
tions in this domain are generally similar to those encountered 
in the frequency domain. Exact subtraction of Xj^(n) from x(n) 
followed by computation of the inverse cepstrum (figure 8 ) 
yields the sequence X2 (n) , exactly, in the time domain. 
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Figure 8 



It is this property which renders homomorphic processing valu- 
able for deconvolution. 

As an example of complex cepstrum transformation consider 
a signal composed of a short pulse (2 samples) convolved v;ith 
a decaying, periodic impulse train. 



4 
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> * 

* 


1 _ 










i 





where 



x(n) = s(n) * p(n) 

cx> 

x(n) = (6(n) + i 6 (n+1) ) * I (-1) (n-kT) 

-k=o 

I R I <1 and T > 2 



Taking the z-transform, 



_ It V -VT 

X(z) = (1 + z/2) ♦ I (-1)' R' z 

k=0 



X(z) = (1 + z/2) • (1 + Rz ^) 



The logarithm then produces a sum 



X(z) = logCl + z/2) - log Cl + Rz ) . 
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Both terms are simply expanded in Laurent series 



log a + z/2) = I z ^ f lz|< 2 



n=“l n2 



„ 00 / n X n^n „ „ 

log (1 + Rz = I z ^ , I Rz ~ | < 1 



n=l 



n 



The Z"transforms 
s(n) = 

p(n) = 

x(n) = 



are easi 
n 



(- 1 ) 



n2 



-n 



(-l)^R^^ 



s (n) + p 



ly recognized. 

, n = -1, ~2 . . .-00 

6 (n~kT) , k = 1,2, 

(n) 



• °° 



V7e see that s (n) occupies only the negative quefrency 
/\ 

region and p(n) only the positive quefrency region. Exact 
deconvolution can be accomplished in this case by zeroing the 
desired half of the cepstrum. 

The canonic form of homomorphic systems for convolution is 
shov/n in figure 9. 




y (n) 



Figure 9 
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The characteristic systems D* and are shown in figures 7 

and 8 respectively. The system L is a conventional linear 
system. When L has a transfer function of 1, then x(n) is 
recovered exactly at the output of The choice of L v/ill 

determine the effectiveness of deconvolution for a given input 
sequence . 

Tv/o further specifications are required to ensure the 
validity and uniqueness of the transformation. The complex 
logarithm is a multivalued function with an infinite number of 
branches . 

log[X(z)] = log|x(z)| + j^Arg[X(z)]± 2Trkj for all k 

v;here Arg specifies the principal value of the phase. This 
ambiguity must be resolved while simultaneously satisfying the 
requirement that X(z) be a valid z-transform. Note that if 
x(n) is to be real and stable, X(z) must be conjugate symmetric 
and analytic in an annulus of the z-plane containing the unit 

/N 

circle. That is, the real and imaginary parts of X(z) must be 
continuous functions of z in the region including the unit 
circle. The imaginary part, arg[X(z)], can only be made 
continuous by "unwrapping" Arg[x(z)] in such a way that all 
jumps of ± 2TTk are removed. This unique unwrapping leads to a 

/S /N 

unique, valid X(z) which transforms to a stable, real x(n). 

The requirement of a continuous phase curve poses some computa- 
tional difficulties which will be discussed in the following 
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section . 

It is apparent from the foregoing description that 
homomorphic systems have the potential for separating convolved 
signals. One might expect, however, by analogy v/ith linear 
systems that deconvolution is most effective for signals with 
certain cepstral properties. This is, in fact, the case and, 
fortunately, seismic signal components are generally amenable 
to deconvolution. Recall that a seismogram is modelled as a 
convolution of a source signature and an impulsive reflector 
series. Reverberation appears as a minimum phase, periodic 
addition to the reflector series. Thus, v;e have 

r^Cn) = p(n) * (b(n) + m.(n)} 

v;here p(n) is the source signature, b(n) is the desired signal 
and m(n) is reverberation. Generalizations can be made concern- 
ing the cepstral properties of each component. 

The source signature is, in general, a mixed phase, short 
duration time sequence. It is clear from the definition of the 
z-transform that any such finite sequence transforms to a 
rational function of z with no poles except at the origin. In 
general 

, m 1 

PCz) = C z n Cl-a.z“^) n tl-b.z) laJ,lb.[< 1 
i=l ^ j=l 

The a^ and b^ represent zeros inside and outside the unit circle 
respectively. z ^ corresponds to a linear phase shift. 
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Here it has been assvuned that the linear phase term is removed 

/s 

before computation. p(n) is a two-sided sequence which is 
always of infinite duration but decays faster than 1/n. 

Hence, most of the cepstral energy is concentrated near the 
quefrency origin. 

The reflector sequence is modelled as a train of randomly 
spaced impulses v/hich may be mixed phase. Stoffa, Buhl and 
Bryan [5] give a general, but very complicated expression 
for the complex cepstrum of such a sequence. Some specific 
examples are given by Schafer [7]. The resulting cepstrum is 
an impulse train v;ith impulses at the time domain impulse 
locations, at all their multiples, and at various other loca- 
tions, both positive and negative on the quefrency axis. Three 
important observations can be made. 

The cepstrum of a minimum phase reflector train contains 
no contributions for negative quefrency. Consider the special 
case of equation (6)- in which all the b^ are zero. This 
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corresponds to all zeros being inside the unit circle, i.e., a 
miniitium phase sequence. Note that p(n) and s(n) in the example 
are minimum and maximum phase respectively, leading to causal 
and anticausal cepstra. 

Secondly, it happens that no non-zero cepstral contributions 
occur betv/een the origin (first impulse) and the location of 
the second impulse in time, if the time series is minimum phase. 
Therefore, the cepstrum of a minimum phase impulse train, unlike 
that of a general sequence does not have its contributions 
concentrated near the origin. The cepstrum of a minimum phase 
impulse train will always contain a gap equal to that between 
the first two time domain contributions. 

Finally, we observe (see Schafer [7]) that a reflector 
series can easily be made minimum phase by exponential v;eighting. 



r' (n) 

R' (z) 



w^r (n) 



Iwl < 1 



C z 



-k 



m /■ w 'i 

n l-(a.w)z l-(b^/w)z 
i=ll 1 J ( 



The value of w is chosen so that 

Ib.w ^1 > 1 for all b.. 

Weighting of the impulse train is effected by w’eighting of the 
entire signal since, if 



s (n) = p (n) * b(n) 
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then 

w^sCn) = (w^p_(n) ) * (v/^b (n) ) . 



Note that the reflector series may be made minimum phase v/ithout 
making the entire signal minimum phase. Very little weighting 
is normally required in practice. 

The remaining contribution to the seismic signal is reverbera- 
tion. This component is merely a special case of a minimum 
phase impulse train in which the impulses are periodic. It is 
easily verified that if 

a 

m(n) = I y(n)6(n-kn ) 
k=l ° 



then 

/S /N 

m(n) = y (n/n^) . 

The cepstrum is, therefore, periodic with the same period as 
the reverberation. 

A useful property for deverberation is derived by Stoffa, 
et al.[5]. The derivation is summarized here because of its 
direct pertinence to seismic processing. 

Consider a normalized multiple signal, 

CO 

ra(n) = I {-l)^K^6(n-i2T ) (7) 

i=0 



where R is the bottom reflection coefficient, 
have seen 



m(n) 



OO 



= I 

i=l 



(-l)^R^ 



6 (n-2iT ) . 
w 



Then^ as we 
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Subtracting the first j cepstral contributions and transforming, 

j ' i 

m. (n) = m(n) - I — <S (n-2iT ) 



w 



i=l 



(z) = I 



„i -2iT 
R z w 



i=j + l 



M. (z) 



= n exp [-- . 

i=j+l i 



Expanding in a power series 



oo oo 



M. (z) = n I - 



„ik -2ikT 
R z w 

TT 



i=j+l k=0 kli 
k-1 



, ^ -2T ,T+X oo 00 -2T 2j+k 

Mj(z) = 1 + I + I I 



j+k 



= 3 i=l ^ (^■'■ 3 ) (j+k-i) 



+ 



(8) 

Comparing (7) and (8), the first j time domain multiples have 
been removed completely and the (j+l)st multiple is reduced by 
l/(j+l) . All succeeding multiples are also reduced. Thus, 
removal of only the first cepstrum multiple would remove the 
first time domain multiple and reduce the second by 1/2, the 
third by 1/3, etc.. 

Having seen the cepstral properties of each seismic signal 
component the advantages of homomorphic deconvolution are 
apparent. The source signature and reflector series have their 
cepstral energy concentrated in different regions of the quefrency 
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axis, thus facilitating removal of the source. The cepstral 
contributions of the reverberation occur at predictable locat- 
tions so that multiple removal is possible. 

Due to the nature of the homomorphic transformation, 
linear filtering of the cepstrum is not the normal convolution 
operation but a simple zeroing of the unwanted contributions. 
That is, the linear filter is generally frequency invariant 
rather than time or quefrency invariant. The name "quefrency" 
was adopted to reflect this reversal of the customary time and 
frequency filtering roles. 

The foregoing analysis is based on assumptions similar to 
those employed in classical seismic processing. Namely, 
we have assumed that the seismogram consists of a source 
signature convolved with distinct, impulsive reflectors and 
periodic multiples. It is difficult to predict the sensitivity 
of the overall processing scheme to these assumptions because 
of its complex analytical structure. Hence, various parameters 
have been varied in the performance analysis to obtain an 
empirical measure of this sensitivity. 

Finally, we note that the additive noise was not included 
in the analytical formulation of the homomorphic processing 
scheme. The algorithm is designed to separate convolved 
components and, unfortunately, no effective processing gain is 
achieved over added noise. In practive, additive noise has 
been dealt with through classical bandpass filtering. This 
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performance analysis includes a description of the effects of 
additive noise on homomorphic deconvolution. 

D . Implementation of Homomorphic Deconvolution 

Figures 7, 8 and 9 show the sequence of operations required 
for homomorphic deconvolution. A processing scheme was designed 
to implement this algorithm in Fortran IV on a 32K digital 
computer. A flow diagram of the scheme is shov/n in figure 10. 

Some aspects of the computation are noteworthy. 

All z-transforms in the algorithm are implemented via FFT. 
Recall from the analytical formulations that z-transforms 
involved in the processing of a real, stable sequence are 
required to have regions of convergence which include the unit 
circle. The discrete Fourier transform is simply a sam.pling 
of the z-transform on the unit circle which, for properly band- 
limited signals, is sufficient to specify the signal completely. 
Data sequences are normally padded with zeros to reduce cepstral 
aliasing, e.g., 2048-point cepstra are computed for 1024-point 
seismograms. Since the cepstrura is alv;ays of infinite extent, 
a truncated version alv/ays results in some aliasing when 
computation is not done recursively. 

The major difficulty in computing an accurate cepstrum is 
the computation of a continuous phase curve. The ta sequences 
are normally sampled at a rate based on the frequency content. 

There is no assurance, hov/ever, that this sampling rate is adequate 

/s 

to uniquely specify X(z) = log[X(z)]. 



SIGNAL x(n) 
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Figure 10 Flow chart of homomorphic dereverberation algorithm. 
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The method used for this analysis is due to Tribolet [9] and 
is thought to be quite accurate and efficient. The flow 
diagram for this algorithm is shov/n in figure 11, 

The phase derivative, and principal value, Arg[X(z)], 

are easily computed from the transforms of x(n) and n x(n) 

(see appendix A) . The values of the derivative are integrated 
using the trapezoidal rule and the integral output is compared 
v/ith Arg X(z) at each step. If the two values do not agree 
within 

2nir ± e n = 0,±1, +2,... 

v/here e is a small positive number, the latest computed value 
of arg X(z), say , is discarded.- The routine then returns 
to the last correct integration value, computes an intermediate 
derivative value, and begins integrating v;ith a step size half 
that of the original grid. The integrate-and-compare process 
is continued at this step size until a^ has been computed 
correctly or until a comparison fails. Integration is resumed 
at the initial step size in the former case or, in the latter, 
step size is again halved. The number of possible step sizes 
is theoretically unlimited. The value of e may be adjusted 
by trial-and-error for most efficient integration. The 
adaptive step feature compensates for the undersampling problem 
in a very efficient and accurate manner. 

Linear phase contributions are easily identified and 
removed from the computed continuous phase curve. Having 
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11 Flow chart of Tribolet phase unwrapping algorithm 
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computed the log magnitude and phase, the inverse z-transform 
is straightforward. 

Linear filtering is accomplished by zeroing the unwanted 
cepstrum values. One might consider windowing procedures which 
are common in linear filtering, but these were not employed in 
the present analysis. 

The inverse cepstrum computation is completely straight- 
forward since no ambiguities arise in the exponentiation 
process. The final steps are shifting the output sequence 
by the linear phase value and unweighting the shifted sequence, 
if necessary. 



- 42 - 



CHAPTER III 

DESCRIPTION OF THE PERFORMANCE ANALYSIS 

The purpose of this analysis is to evaluate, both quantita- 
tively and qualitatively, the performance of the TDL and homo- 
morphic dereverberation techniques. Each method is to be 
evaluated for a variety of simulated seismic conditions in order 
to determine those factors which significantly influence perfor- 
mance. The comparative nature of the analysis is intended to 
emphasize the relative strengths and weaknesses of each technique. 

It should be noted that absolute performance figures are 
not inherently valuable, especially when obtained from synthetic 
data. The diverse geological and oceanographic conditions 
encountered in marine seismology coupled with the many different 
processing systems currently employed may be expected to yield 
a range of absolute results. The greater value of this analysis 
is to indicate the parameters, environmental and mathematical, 
which can be expected to affect significantly the performance 
of these algorithms. The numbers obtained provide a measure 
of the relative performance of the two methods under similar 
conditions and, in some cases, provide asymptotic performance 
criteria. Synthetic data were chosen so that signal parameters 
could be accurately controlled. 

Three criteria are specified for comprehensive evaluation 
of performance. The first, most direct, measure of effectiveness 
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is the percentage of multiple energy removed from the signal. 

This is easily measured by calculating the zero-lag correlation 
of the signal, in time windows spanning only the multiple 
locations, before and after processing. That is, the squared 
amplitude (energy) of the signal in the multiple region is 
computed after processing and subtracted from the squared 
amplitude of the original signal in that region. This 
difference is divided by the original energy of the multiple 
to yield the fraction of energy removed by processing. The 
use of synthetic data facilitates this method of analysis 
since reflectors and multiples can be placed in disjoint regions 
to avoid ambiguity. 

The second criterion is the amount of signal distortion 
introduced by dereverberation. This is measured by com.paring 
reflector energy before and after processing. As before, 
reflectors and multiples must be disjoint for meaningful results. 

Although these tv;o criteria provide an accurate quantitative 
measure of performance they are restricted to situations in 
v;hich reflectors and multiples do not overlap. The overlap 
case is most important in processing real data since the 
multiples then abscure the reflectors most severely. In order 
to judge performance in these situations we must evaluate 
qualitatively the improvement of visual interpretability . 

This visual enhancement of ref lector-to-multiple ratio is our 
third criterion. It is an important measure in spite of its 
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subjective nature since the primary means of seismogram 
analysis is visual interpretation. 

The analysis is limited in scope to the single channel 
processing configuration. Although both techniques are 
applicable, in principle, to multichannel processing, the many 
additional variables involved and unavailability of appropriate 
synthetic data would lead to a complicated extension of this 
analysis. A single channel treatment is adequate to identify 
the important performance traits of both methods. 

Parameters to be varied fall into two general categories; 
environmental and operational. The environomental parameters 
include noise level, multiple periodicity, multiple-to-signal 
level and multiple/ref lector separation or overlap. These are 
varied within ranges v/hich are . thought -to be representative of 
ambient conditions normally encountered in marine seismology. 
Effects of noise level are considered only for the case of 
white Gaussian noise. The effects of aperiodicity have not 
been evaluated for the TDL algorithm because it is intended 
primarily for deep v/ater use where only the first multiple is 
usually of interest. Operational parameters refer to those 
V7hich can be controlled during processing. These .include filter 
cutoff frequencies, operator lengths, travel time estimates, 
cepstral stopbands, and v;eighting. Windowing of the correlation 
function is not evaluated although a discussion of this subject 
is contained in reference [3J . 
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The algorithm used in generating simulated earth impulse 
responses is due to Theriault [lOJ . A brief description of 
Theriault's earth model is given here. 

The model is based on the following assumptions : 

Cl) An air gun source generates longitudinal pressure 
waves w’hich impinge upon all earth layers as 
normally incident plane waves. 

(2) All earth layers are horizontally infinite, 
parallel and homogeneous. 

(3) Abrupt changes in acoustic impedance occur at 
each layer interface and these boundaries are 
characterized by the customary acoustic reflec- 
ion and transmission coefficients. 

(4) The earth has a uniform density. 

(5) The water column is a non-attenuating fluid 
v/ith a perfect pressure release interface at 
its surface. 

(6) Each layer is characterized by a transfer 
function, F(w) v;hich represents the attenuation 
and travel time delay for that layer. 

These assuraptions are incorporated into a lujnped parameter 
model. Figure 12 shows a frequency domain model of a two 
layer earth. The Fj^^Cw) have the functional form 
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Figure 12 Two layer earth transfer function schematic 
(after Theriault) . 
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where x^, and are the ith layer thickness, sound speed 
and layer attenuation parameter respectively. These transfer 
functions are combined using a semi-group property rather than 
the usual frequency domain multiplication. 

The -1 multiplier completes a feedback loop around the 
source which generates the water column multiples. Since the 
water column is assumed to have no attenuation the multiples 
appear in the earth response as impulses, and in the resulting 
seismogram as replicas of the source signature reduced by the 
bottom reflection coefficient. 

The above multiple mechanism is inadequate for representing 
effects of incoherent (distorted) multiples and varying bottom 
interaction mechanisms. These effects are introduced by inser- 
tion of v/ater column attenuation v;hich causes the bottom response 
and multiple responses to be of exponential form given by the 
Fourier transform of (9) . Thus the bottom response is extended 
in time and each multiple is a distorted version of the previous 
one. Examples of multiple appearance v:ith and without water 
column attenuation are shown in figure 13a and b. 

The topmost loop of figure 12 simulates the effects of 
finite receiver depth x. Any number of layers with the desired 
parameters can be combined into an overall earth transfer 
function which is easily transformed to yield the earth impulse 
response . 

Synthetic seismograms for this analysis were generated by 



(a) 







Figure 13 (a) Synthetic Pei^r'.ogran v/ith coherent multiples. 

(b) Synthetic seismogram vaith distorted multiples 
due to water column attenuation. 
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convolution of various Theriault impulse responses v;ith air 
gun signatures obtained from actual at-sea recordings. A 
typical signature is shov/n in figure 14. 

It is important to note that these seismograms contain 
all water column multiples and internal multiples of all 
orders. For realistic parameter values earth attenuation 
characteristics usually render internal m.ultiples negligible. 

Finally, v;e note one drawback of using this model for 
multiple removal analysis. The algorithm produces multiples 
v;hich are exactly periodic. This periodicity gives these 
seismogram.s strong correlation characteristics which are not 
usually encountered in practice. The lack of periodicity in 
actual seismograms is due to the horizontal separation of the 
source and receiving array. The difference in travel paths 
arising from this separation is illustrated below for a primary 
return and first multiple. 




multiple path length = 



a > 8 
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Figure 14 Air gun signature (obtained by at-sea recordings) 
convolved with synthetic earth response functions 
to form seismograms. 



- 51 - 



For a knovm separation and water depth the travel time 
difference can be easily ca’lculated. This effect becomes 
minimal in deep water where cos a and cos 3 are approximately 
equal to one . 
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CHAPTER IV 
RESULTS 



A. Introduction 

The results of applying the multiple removal methods 
described in Chapter II to synthetic data are described here in 
terms of the criteria of Chapter III. Performance based on the 
first two criteria is emphasized because it can be quantified 
much more accurately. Specifically, the reflectors and 
multiples have been positioned at distinct locations in most 
cases so that the amount of energy removed from reflectors and 
m.ultiples can be measured without ambiguity. This emphasis 
leads, hov/ever, to a certain lack of realism in several of the 
synthetic data plots. Separation of this kind in an actual 
seismogram would, of course, eliminate the necessity for 
multiple removal. Hence, several cases of interfering reflectors 
and multiples are also shown. Although these are not amenable 
to quantitative analysis they can be judged on the basis of 
the third criterion, viz., improvement of visual record quality. 

A second deviation from normal processing conditions has 
been required to compare effectively the performance of the 
two methods. Homomorphic dereverberation is accomplished in 
practice Csee {5J ) concurrently with source deconvolution, 
where practical, since both operations are simply performed 
after the cepstrum has been computed. This leads to a 
considerable amount of energy removal at each multiple and 
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reflector location which is due to source deconvolution alone. 
Although this may lead to improvement of the record, the criterion 
of energy removed does not accurately measure de reverberation 
performance in the same way it does for the TDL results. For 
example, source deconvolution might typically lead to 90% 
reduction of multiples and 90% removal of reflector energy. 

This gross change in configuration cannot be effectively 
compared with TDL processing on the basis of energy removed. 

Hence, source removal has not been accomplished in most cases 
tested. The results of this "partial" processing can be 
quantitatively evaluated and easily compared with TDL results. 
Cepstral filtering v/hich includes source removal, thus retain- 
ing only high quefrency energy, is referred to as "longpass 
filtering". Several examples of longpass filtering are included 
and interpreted in terms of the third criterion. 

The performance of each method is discussed individually 
for various processing situations. A direct comparison of the 
performance of both methods for the same data is also included. 

The comparison is extended in Chapter V. 

B . Results of TDL Dereverberation Performance 

1. Operator Length, Multiple Distortion and Multiple-to- 
Signal Ratio 

The number of tap gains required for optimum multiple 
removal was found to be highly signal-dependent. Recall from 
Chapter II that the tapped delay line is essentially an 
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estimate of the high energy portion of hj^(t), the earth impulse 
response. In most applications the reflected bottom response 
is dominant. The time-bandwidth product of this response 
would then be expected to govern the filter length requirements. 
Measured results confirm this. 

Seismograms containing exactly impulsive (one sample) 
bottom responses exhibited little or no variation of performance 
with filter length. Typical filter impulse responses for a 
seismogram of this type are shown in figure 15. The reflection 
coefficient in this case is 0.3. The first tap gain is close to 
-0.3 in each response, as would be expected from the Backus 
formulation in which the second operator point is an estimate 
of the bottom reflection coefficient. Increasing filter length 
can be seen to cause variation in the "estimation" of the reflec- 
tion coefficient. Figure 16 shows energy removed vs. filter 
length for this seismogram. Multiple energy removed decreases 
slightly with increasing filter length. Reflector distortion 
is nearly constant at lov/ values (6% for one and -2% for the 
other). The signal used in figure 16 is shown in figure 17a. 

The processed result shov/n was obtained using only one tap. 

Introduction of a non-impulsive multiple mechanism was 
found to produce a marked dependence of performance on operator 
length. Synthetic seismograms were generated with finite 
length bottom responses, resulting in distorted multiples 
which are not simply weighted replicas of the source signature. 
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Figure 15 TDL impulse responses for three operator lengths, 
(a) 5 taps (b) 10 taps (c) 15 taps 
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Figure 16 Operator length vs. perfomance for a signal with 
v/ith impulsive multiples. 
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Figure 17 (a) Seisriograin with coherent i^ultiples. 

(b) Result of TDL processing with a one point 
operator. 



-58- 



Figure 18 shows filter perforraance vs. operator length 
for three seismograms having different degrees of distortion 
in their bottom reflection mechanism.s. This distortion is 
equivalent to extension of the bottom impulse response in time. 
The signals associated with curves (1) , C2) and C3) are shov/n 

in figures 19a, 19b and 19c respectively. The increasing 
distortion of the multiple at 2.0 seconds is evident. Curve (1) 
corresponds to a signal with very slight multiple distortion. 
Operator length is seen to have no effect on performance. 

The seismogram corresponding to curve (2) has a bottom impulse 
response which is significant for T = .03 seconds. The tap 
spacing in this case is 1/2W = 4.88 msec, so that 2TW = 6.1, 
and six or seven tap gains should be adequate if the signal 
has been properly sampled. Reference to curve (2) confirms 
that increasing the filter length beyond seven does not improve 
performance. Filters of fewer than seven elements yield 
monotonically decreasing performance. The bottom response 
for curve (3) is significant for .045 seconds so that, for the 
same tap spacing, 2TW = 9.2 and we anticipate that nine or ten 
taps v;ill be adequate. This is, in fact,, the case. 

Figure 20 shows the 5, 9 and 15 point filter impulse 
responses associated with curve C3) . There is an observable 
convergence to a t e ^ shape which is the actual functional 
form of the synthetic bottom response. Figures 20a and 20b 
exhibit a "diverging tail” effect which was found to be common 
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Figure 18 Performance vs . operator length for three signals 
with bottom interaction times varying from (1) 
impulsive to (3) .045 seconds. 
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Figure 19 Seisinograns associated with the curves of figure 

18. (a) Curve (1). (b) Curve (2). (c) Curve (3). 
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Figure 20 Filter impulse re?^ponses associated v'ith Curve (3) 
of figure 18. (a) 5 points (b) 9 points 

(c) 15 points. 



-62- 



when the specified filter length is too short. This type of 
phenomenon occurs in some numerical approximation methods when 
an inadequate nximber of terms is specified. 

It is apparent from the above results that the effect of 
operator length is closely related to the bottom interaction 
mechanism. As discussed in Chapter II, the filter should be 
an estim.ated replica of the bottom impulse response when the 
interaction process can be accurately modelled as a convolution 
of the source signature with the bottom response. Figures 
15 and 20 are good excimples of this behavior. 

In the cases summarized in figure 18 performance increases 
as multiple distortion increases. This need not be true in 

general since bottom interactions nay become very complex. 

2 ""b t 

The slov/ly varying f^e responses lead to operators V7hich 
have a greater cancellation effect as they are extended. A 
higher bandwidth bottom impulse response might not exhibit 
this behavior. As it was not possible to include more complicated 
bottom responses in the earth model used, these effects were 
not investigated further. 

Reflector distortion was found to be nearly constant for 
all filter lengths tested. The reflector at 2 . 7 seconds v;as 
essentially undistorted in all cases. The 3.5 second reflector 
had an average distortion of 7%. Figures 21-24 show some 
processed results for the signals in figure 19 . Each of the 
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Figure 21 (a) Seismogram with very coherent multiples. 

(b) Result of TDL processing with three taps. 
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Figure 22 (a) Seismogram with moderately distorted multiples; 

bottom response is significant for .035 seconds, 
(b) Result of processing v/ith 7 taps. 
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(a) Seis"^ogram v/ith considerable multiple distortion; 
bottom response is significant for .045 seconds. 

(b) Result of TDL processing with 11 taps. 



Figure 23 
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first three figures C21, 22 and 23) show signals before and 
after processing with the fev/est number of taps required to 
achieve optimxun performance. Visually, multiple removal is 
almost complete in each case. The signal of figure 23a is 
shov/n in figure 24 after processing v/ith 3, 5 and 7 taps. 

The improvement from figure 24a to 24c is obvious. 

Multiple-to-signal ratio (defined here as the ratio of 
energy in the first multiple to that in the largest sub-bottom 
reflector, abbreviated MSR) was found to be of little importance 
in most cases of interest. For signals in which multiples are 
large enough to be a problem (comparable to, or larger than 
smaller reflectors) , the multiple dominate the crosscorrelation 
function so that an effective filter is generated. For these 
cases the performance v’as found to be insensitive to Uie v/idth 
of the time window used for correlation. In the relatively 
less interesting case of signals with small multiples, the 
performance is greatly dependent on the choice of correlation 
window. If large reflectors are included in the window, perfor- 
mance is adversely affected because of the large reflector 
contribution to the statistics. For some cases of interest 
this effect may be a consideration in choosing an appropriate 
correlation window. Inclusion of large reflectors should be 



avoided . 
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Figure 24 Results of processing seismogram of figure 23a v;ith 

inadequate TDL lengths (a) 3 taps (b) 5 taps (c) 7 taps. 
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2. Water Travel Tirae Estimate 

It was noted in Chapter I.J .that an estimate of the water 

colximn travel time is required for implementation of the TDL 

algorithm. This is accomplished in practice by various methods 

including visual estimation, energy detection and correlation 

techniques. The estimate appears in the filter design equations 

as the minimum shift of the crosscorrelation function, R 

mr 

This estimate may also be identified as the prediction distance 
when the operator is interpreted as a prediction error filter. 

The actual estimation of water column travel time was not 
investigated in this analysis. The effects of travel time 
estimation on filter performance were, however, considered. 

Figure 25 shov.'s the results of water travel time estimation 
errors for the ideal case of a signal with impulsive multiples. 
The unprocessed seismogram, shown in figure 26a, contains 
reflectors at 2 . 7 and 3.5 seconds and a strong multiple at 2.0 
seconds. The actual tv;o-way travel time in this case is 1.0 
second. The strong similarity between the bottom reflection 
and multiple is apparent. 

First multiple energy removed is very sensitive to travel 
time estim.ation, v/ith an error of ± 10 msec resulting in a 
performance degradation of about 20% . This effect is analogous 
to that observed in matched filter receivers in that the 
coherent signal components exhibit very high correlation for a 
small range of lags.' In this case the strong coherence and 
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Effects of V7ater coluran travel tiine estimate on 
multiple removed for a signal with coherent multiples. 



Figure 25 
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Figure 2 6 (a) Unproces<^ed seismogram v;ith impulsive multiples. 

(b) TDL result based on correct water travel time 
estimation . 
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Figure 26 (coni' d) (c) Result of ']?DL processing based on an 

early water travel tirne estinate. 

(c) Same processing based on a late estimate. 
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raultiple periodicity yield the sharp correlation peak for lags 
near the water colvunn travel time, which is also the multiple 
period . 

The second and third multiple performance peaks occur 
at 1.08 and 1.04 seconds respectively, and have values of 24% 
and 60%. These shifted and reduced peaks are caused by 
digitization effects. For such an extremely coherent signal, 
a multiple position deviation of one sample Cdue to sampling 
interval round-off) can cause the cancellation operation 
to be severely affected. In this case the second and third 
multiple performance maxima are due to secondary crosscorrelation 
peaks introduced by the periodic oscillations in the source 
signature. These peaks occur at intervals approximately 
equal to the period of the basic source (air gun) frequency. 

The three performance peak values correspond closely to the first 
three air gun pulse amplitudes. 

Figures 26 b, c and d show the results of on-time, early 
and late travel time estimates respectively. In the first case 
the multiple has been effectively removed while the early and 
late estimates lead to multiples which are still significant 
after processing. 

The effects of travel time estimation error on the 
operator impulse response are seen in figure 27. The optimum 
estim.ate yields an operator with a large peak at the origin, 
nearly equal to the bottom reflection coefficient (.2) , and a 
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(a) 





Figure 27 TDL filter impulse responses based on travel time 

estimates which are (a) 40 msec early, (b) correct, 
and (c) 20 msec late. 
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smaller peak at the source pulse oscillation interval. The 
initial spike is due to the largest peak in the crosscorrelation 
function Capproxiraated here by the correlation function) v/hich 
occurs at a shift equal to the two-way travel time. The smaller 
peak in the filter is introduced by an additional shift of one 
source period. The early (by 40 msec) estimate still allows 
observation of the crosscorrelation maximum and yields the 
shifted spike of figure 27a. It is evident in figure 25 that 
performance is reasonably good for estimation errors less than 
about .07 seconds, which is the length of the tapped delay 
line. For greater travel time errors the range of correlation 
lags' computed does not include the peak; hence, -the filter is 
ineffective. Late travel time always results in poor perfor- 
mance since the peak is not observed. Figure 27c shov;s a 
typical impulse response due to late travel time estimation. 

The seismogram of figure 28a contains a small reflector 
at 2,2 seconds and a larger one at 3.0 seconds. The first 
multiple partially overlaps the smaller reflector and the second 
multiple coincides v/ith the larger. Figures 28 b, c and d 
demonstrate how travel time estimation errors can lead to 
ambiguous results. Figure 28b was processed using the correct 
travel time, resulting in good resolution of both reflectors. 

The early and late estimates lead to the signals of figures 
28c and d in which the smaller reflector cannot be clearly 



resolved . 
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(a) Seismogram with overlapping reflectors and 
multiples . 

(b) Result of TDL processing based on a correct 
travel time estimate. 



Figure 28 
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Figure 28 (cont'd) (c) Result of TDL processing based on a 

40 msec early travel time estimate. 

(d) Same processing based on a 20 msec late 
estimate . 
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Less coherent seismograms exhibit less sensitivity to 
travel time estimation as shown in figure 29, Curves (1) , 

C2) and (3) were generated using the seismograms of figure 19 a, 
b and c respectively, which contain increasing distortion in 
their multiples. The MSR is considerably lower than in 
figure 25 so that the later multiples are not clearly visible. 

The crosscorrelation peaks are considerably broader than in 
the impulsive multiple signals. Curve (3), the least coherent 
signal, exhibits the lowest sensitivity to travel time estimates. 
Increasingly sharper peaks are evident for curves (2) and (1) . 

The performance degradation due to over-estimation is precipitous 
in all cases although the rate of fall-off is related to 
signal coherence . 

Reflector distortion (not shown) was negligible in all 
cases for the first reflector and nearly constant at 7% for 
the second . 

Examples of processed signals for curves (2) and (3) 
are shov.Ti in figures 30 and 31. 

3. Multiple Periodicity 

The effects of aperiodic multiples on TDL dereverberation 
performance were not investigated in this analysis since the 
algorithm is designed primarily for deep water use, where only 
one multiple is significant in most applications. Successful 
employment of this technique in shallow water is limited since 
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Figure 29 Effects of water column travel time on three signals 
with varying degrees of multiple distortion. 
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Figure 30 Processed oeisnograris associated with figure 19b 
and curve (2) of figure 29. (a) Travel time 

estimate 40 msec early (b) 40 msec late. 
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Figure 31 Processed scisr.ograms associated with figure 19c and 
curve (3) of figure 29. (a) Travel tine estiinate 40 

.msec early (b) 40 msec late. 
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the approximation of the crosscorrelation, by the 

correlation, near the multiple onset is not generally 

valid. Reflector energy in the multiple regions is usually 
significant in shallow water seismograms. 

It is apparent from the TDL formulation and the discussions 
of the preceding tv;o sections that the TDL algorithm has 
some capability to compensate for aperiodicity in which the 
deviation is .of the order of the tapped delay line length. 

For greater aperiodic! ties the TDL filter v;ill not be effective 
on later multiples. Figure 25 indicates that even slight 
deviations can be very detrimental in later multiple removal. 
Dynamic corrections can be applied as suggested by Backus [1] , 
however, these are beyond the scope of this study. 

4. Reflector/Multiple Overlap 

Figures 32-36 show seismograms before and after processing 
for several cases in which multiples and reflectors overlap. 

A brief description of each situation is given here. 

The seismogram of figure 32a has reflectors at 1.5, 2.1 
and 3.0 seconds, and distorted multiples. Significant energy 
near 2.0 seconds is removed by processing but a clearly visible 
response is still present. Other regions of the signal are 
not visibly affected. The visual resolution is not significantly 
improved in this case, however, stacking of successive shots 
after multiple removal could be employed to enhance the reflector 
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Figure 32 Seismogra-m with reflectors at 1.5, 2.1 and 3.0 

seconds (a) before and (b) after TDL processing. 
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in this situation. The success of stacking would require 
that tile reflector be more coherent than the residual multiple 
after dereverberation. 

In figure 33 reflectors occur at 1.5, 2.02 and 3.0 seconds. 
Reflectors may be expected to have a strong influence on the 
crosscorrelation function in this situation of extreme overlap. 
The result is nearly complete removal of the second reflector 
with the multiple. The last reflector is clearly visible at 
3.0 seconds. In this example the tapped delay line length 
C.054 second) is greater than the reflector/multiple separation 
C.02 second). The bottom impulse response is significant for 
.035 second so that it cannot be estimated accurately with an 
operator which is shorter than the reflector/multiple separation. 
Even if the bottom response v^ere much shorter, the effect of 
the large reflector in the crosscorrelation windov/ would lead 
to degraded performance. This example shows worse degradation 
than would be expected in practice because the reflectors are 
very coherent in this synthetic data. The less coherent 
reflections in deep v/ater signals v;ould generally be less 
susceptable to removal. When reflector energy is significant 
in the crosscorrelation, hov'ever, and TDL length is greater than 
reflector/multiple separation, the filter design algorithm 
essentially interprets the reflector as part of the multiple 
and tries to remove it. 

Figure 34 contains reflectors at 1.6, 1.95 and 3.0 seconds. 
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Figure 33 Seisnogran v;ith reflectors at 1.5, 2.02 and 3.0 
seconds (a) before and (b) after TDL processing. 
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Figure 34 Seismogram with reflectors at 1.6, 1.95 and 3.0 
seconds (a) before and (b) after TDL processing. 
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These locations a.re slightly different from those of figure 33. 
Significant energy is present at the second reflector location 
after processing. Since the reflector precedes the first 
multiple in this case it has a less significant effect on the 
crosscorrelatioa f;anction than the 2.02 second reflector in the 
previous example. The variability of performance on individual 
seismograms (of similar structure) with overlap again suggests 
the use of stacking to enhance reflectors. 

The reflectors of figure 35 are situated as in figure 34; 
however, the first reflector and multiple are significantly larger 
in this case. This seismogram is more typical of a shallow 
water return. The processing is quite effective on the first 
and second multiples. The aperiodicity of reverberation in 
most shallow water seismograms would lead to much poorer 
performance on later multiples, in practice. 

Figure 36 illustrates a situation where the overlap is 
not severe but visual resolution is impaired by the first and 
second multiples. Reflectors occur at 1.55, 2.4 and 3.25 
seconds. Processing results in effective reduction of both 
multiples and significantly improved resolution in the second 
and third reflectors. 

We conclude -from these examples of ref lector/raultiple 
overlap that visual improvement due to dereverberation varies 
widely, depending upon several aspects of the individual signal 
structure. The bottom impulse response, the extent of reflector/ 
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Figure 35 Seismogram with reflectors at 1.6, 1.95 and 3.0 
seconds (a) before and (b) after TDL processing. 
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Figure 36 Seismogram with reflectors at 1.55, 2.4 and 3.25 
seconds (a) before and (b) after TDL processing. 
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multiple separation, the relative coherence of reflector and 
multiples, and tfie relative sizes of signal components all 
appear to affect performance. Each of these factors influences 
the effectiveness of the crosscorrelation operation in identify- 
ing energy which is specifically due to the multiples. 

The visual Improvement due to de reverberation can best be 
determined by analysis of continuous seismic profiles (presenta- 
tions showing man.y seismograms side-by-side) since coherent 
residual energy, if present, becomes apparent as visual trends 
in the data. The effects of stacking can also be observed in 
this format. It was not practical to generate such a profile 
with synthetic data but the single-shot results indicate that 
TDL dereverberation is potentially useful for enhancing the 
visibility of reflectors which are partially masked by 
multiples . 

5. Additive IThite Noise and Filtering 

White Gaussian noise was generated in the following manner. 
First, a set of uniformly distributed random numbers was 
obtained using a standard digital routine. A set of approximately 
Gaussian numbers was then formed by summing separate groups of 
twelve of the uniformly distributed numbers. The resulting set 
was weighted to obtain the desired standard deviation. Figure 
37 show’s a seismogram v;ith tv/o different levels of added noise, 
each lowpass filtered at 50 Hz. 
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Figure 37 Seismograms with added white noise, lowpass at 50 Hz. 

(a) a =50 (b) a =100. 
n n 
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Figure 38 illustrates the performance of the TDL filter 
for various noise levels on a seismogram having disjoint 
reflectors and multiples. The maximum noise standard deviation 
shov/n corresponds to a SNR of 0.7 dB. No bandpass filtering 
was done, prior to dereverberation. The results presented here 
are typical of the cases tested. Multiple energy removed 
decreases monotonically v/ith increasing noise level; almost 
linearly in the case of the first multiple. The later multiples 
are more severely affected. As the noise level becomes comparable 
to their am.plitudes, the filter becomes ineffective in removing 
them. This behavior is due to the effect of increasing incoherent 
energy in the signal. The correlation-cancellation operation 
is designed to detect coherent, periodic components. These 
become increasingly masked as more noise is added. For this 
reason seismograms with larger multiples exhibit less sensitivity 
to noise. Seismograms of similar structure (i.e., exactly the 
same reflector and multiple locations) whose multiple energies 
differed by a factor of tv;enty were found to exhibit consider- 
able dereverberation performance differences in high noise. 
Processing of the signal with the larger .multiples resulted in 
25% greater removal of first multiple energy. 

Lowpass filtering before de reverberation can lead to a 
significant improvement in TDL filter performance in some 
signals. The plotted results of lov.'pass filtering two noisy 
signals at various frequencies are shown in figure 39. 
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STANDARD DEVIATION OF ADDED WHITE NOISE 



Figure 3 8 



Effects of noise level on TDL performance for a 
signal with coherent periodic multiples. 
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FILTER CUTOFF FREQUENCY (Hz) 



Figure 39 TDL performance vs. lov;pass filter cutoff frequency 

for signals with coherent (lower curve) and distorted 
(upper curve) multiples. 



-94- 



(A third order Butterv/orth filter was used throughout. The 
visual difference in two filtered signals is shown in figure 40. 
The first is lov/pass filtered at 90 H 2 and the second at 30 Hz.) 
The seismogram used to generate the upper curve contains 
multiples vmich are considerably distorted while the lower 
curve is due to a more coherent signal. Decreasing the pass- 
band from 90 to 30 Hz produces a 12% improvement ,in performance 
in the second case but the less coherent signal is relatively 
insensitive to filtering. Several other signals evaluated 
were found to exhibit the same behavior. The phenomenon is 
apparently due to the averaging, or smoothing effect of the 
operator in the distortion case. Recall that signals of this 
kind tend to have more extended waveforms in their TDL operators 
(See figure 20). Convolution of such functions with a noisy 
signal "smoothes out" the noise because of the significant 
operator length and the random fluctuation of the noise. This 
has the tv;ofold result of better overall performance in noise 
and less sensitivity to removal of higher frequency noise energy. 

Performance degenerates for filter cutoff frequencies below 
30 Hz because the spectral energy of the. signal itself is 
concentrated in this range. Figure 41 shows the distribution 
of energy in the frequency domain for the seismogram associated 
with the upper curve of figure 39. Most of the energy is 
concentrated between 5 and 30 Hz. 

Figures 42-45 illustrate the visual improvement of noisy, 
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Figure 40 Seismograms with equal white noise levels, lowpass 
filtered at (a) 90 Hz and (b) 30 Hz. 
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Figure 41 Squared magnitude of the Fourier transform for the 

signal associated with the upper curve in figure 39. 
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Figure 42 



(a) Unprocessed seismogram with added noise, F^-50 Hz. 

(b) Result of TDL processing. 
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Figure 43 (a) Seismogram of figure 42a with high noise, F^=50 Hz. 
(b) Result of TDL processing. 
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Figure 44 (a) Seismogram with low noise, F =50 Hz. 

(b) Result of TDL processing. ^ 
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Figure 45 (a) Seismogram of figure 44 with high noise, F =50 Hz. 

(b) Result of TDL processing. 
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Figure 46 Operator impulse responses for three different 

noise levels on same signal with coherent multiples. 
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lowpass filtered C50Hz) signals after processing. Two different 
seismograms are shown, each at two different noise levels, 
before and after processing. It can be seen from these figures 
that no serious reflector distortion occurs due to noise or 
filtering. 

The presence of white noise has a predictable effect on 
the appearance of the TDL operator. The impulse response itself 
becomes noisy due to the added incoherent energy and a bias is 
introduced into the large cancellation peak. An example of this 
is shown in figure 46 for a signal with coherent multiples. 

The variation in the first tap gain is significant and the change 
in appearance of the overall impulse response is appearent. 

C. Results of Homomorphic Dereverberation Performance 

1. Introduction 

The theory of homomorphic de reverberation is based on the 
properties of periodic minimum phase impulse train cepstra. The 
incorporation of more realistic effects, such as aperiodicity , 
distorted multiples, and non-impulsive reflectors, leads to 
analytical intractability in most cases. Hence, in presenting 
the experimental results, we view the more complex signal 
structures in terms of their relation to the ideal signal models 
of Chapter II. In so doing we hope to provide a basic reference 
for interpreting observed phenomena, and to provide a stronger 
intuitive picture of the mechanisms at work in homomorphic 



dereverberation . 
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All cepstra used in this analysis are 2048 points in 
length. Seismograms contain 1024 points except in cases where 
resampling was required. Zeros have been appended as necessary. 
The 2048 point length has been used throughout the analysis for 
uniformity and reduction of aliasing. In most practical seismic 
processing, shorter cepstrum lengths are adequate. 

2. Multiple Periodicity 

It was shown in Chapter II that later components of a 
periodic, minimum phase impulse train can be significantly 
reduced by removal of the first one or two cepstral contributions. 
This property is potentially valuable for dereverberation, 
especially in shallov; water cases where several strong multiples 
may appear in the data. In order for this property to be useful 
it must be relatively insensitive to (at least) slight 
aperiodicity since actual reverberation is not exactly periodic. 
The result of Stoffa, et al. summarized in Chapter II is 
extended here to rapidly decaying, aperiodic, minim.um phase 
impulse trains. 

Consider removing the first j cepstral contributions of 



mCn) , the cepstrum of a minimum phase impulse train with 
contributions (,-R) ^ at n^ , i = 1,2,3,..., and |r| < 1. 

We then have 



/ • • • / 






i=j+l 



6 (n-n^) . 
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Trans forming. 



/S / T-»\ ^ 

M.(z) = I - icSL 2-"]i . 
Jl=j + 1 



Exponentiating , 



Mj Cz) = exp 



(X) 0 -- 

^ % 

il=j+l ^ 



M. Cz) = n exp 
^ £=j+l 






This may be expressed as a power series 



M (z) = n ^ .. 

^ £=j + l i=0 i!Jl^ 



00 Q\ 



M.(z) - n I (-R)*^ „-n^ ^ (-R)2*- _-2n, (-R) ” _-3n 

3 t=j+lL" 1 " 2^2 " " 

( 10 ) 

The rapidly decaying coefficients of z combine multiplica- 

tively to yield the time domain impulse coefficients of m^ (n) . 
If only the first cepstral impulse is removed (j=l) , the 
largest value attained by the third coefficient in brackets is 




C-R) 

21^ 



21 



= ^ 
A=2 8 



which is insignificant for reflection coefficients of interest. 
All succeeding terms in CIO) will be vanishingly small. We 
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then have the approximate expression 
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Expanding this product v/e obtain 



= 1 + I 

^ k=l 
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which transforms to 

m . (n) = 6 (n) + • 6 (n-n . , ) + — - • 6 (n-n . 

J .X-. 1+1 j + 2 3^2 



j + 1 



(11) 



It is apparent that the remaining terms have been reduced by 
factors equal to those obtained in the periodic case. 

The above result has been verified experimentally by process 
ing a periodic signal containing only a strong bottom reflection 
(R=.55) and multiples (figure 47a). Removal of the first 
multiple component in the cepstrum results in 75% and 85% energy 
reductions of the second and third multiples respectively. 

These correspond exactly to the 50% and 67% amplitude reductions 
predicted by (11) • The periodicity was then disturbed by 
average deviations of 1%, 5%, 10% and 20% of the original 
period. In each case the later multiple reductions coincided 
with (11) . Figure 47 shows waveforms before and after process- 
ing for the exactly periodic signal and the case of 20% 
aperiodicity . 
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(a) Signal composed of only periodic multiples of 
the bottom response. 

(b) Result of removing the first multiple cepstrum 
contribution. 



Figure 47 
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■ Figure 47 (cont'd) (c) Signal of (a) with periodicity disturbed 

by average deviations of 20% (.2 sec). 

(d) Result of removing the first multiple 
cepstrum contribution. 



-108- 



3. Multiple-to-Signal Ratio 

The results obtained for aperiodic reverberation must be 
extended to signals containing reflectors as well as multiples 
if they are to be of practical use. This is not feasible 
analytically due to the nonlinearity of the homomorphic trans- 
formation and the reflector characteristics in the earth model 

used. Recall that, for this analysis, reflectors have the 
2 “"]d t 

form t e rather than simply impulses. Each seismogram con- 
tains an additive combination of these reflectors with 
multiples. The logarithmic operation on the Fourier transform 
of this sum causes the cepstral contributions due to reflectors 
and multiples to be analytically indistinguishable. Hence, 
this phenomenon has been investigated empirically for various 
reflector sizes. It v;as found that cepstral properties of 
multiples are essentially preserved in the presence of non- 
impulsive reflectors although important effects were evident 
in the cases tested. 

Percentage removal of the first multiple was found to be 
dependent on the width of the cepstral stopband. Figure 48 
Illustrates this effect for seismograms v;ith different reflector 
stengths and periodic multiples. Each signal requires a 
stopband of about 200 msec (41 samples) for complete removal 
of the first multiple. As notch width is decreased the perfor- 
mance becomes increasingly sensitive to MSR. For the smallest 
notch width shown (40 msec) , performance varies almost 30% 
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Figure 48 Effect of cepstral stopband v;idth on homomorphic 

dereverberation performance for signals of varying 
MSR. 
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with MSR. In the absence of reflectors the first multiple 
could be completely removed by zeroing one sample of the cepstrum. 
It appears that the inclusion of reflectors has the effect of 
distributing the multiple energy in the quefrency domain, 
although, in all cases tested the energy remains concentrated 
near the time domain multiple location. Figure 49 shows the 
0.5 second section of the cepstrum centered at the multiple 
location for each of the seismograms of figure 48. The cepstra 
are identical in form but the absolute cepstral energy increases 
with MSR. Three large peaks occur in each cepstriam between 
.94 and 1.0 second. This similarity in form suggests that 
equal stopbands should produce equivalent percentage reduction 
of multiples. We conclude, however, from the experimental 
results that a greater portion of the multiple energy is 
concentrated near 1.0 second in the higher MSR cepstra. This 
effect makes it difficult to select stopband limits by peak 
detection or visual inspection of the cepstrum. Notch widths 
which yield the best trade-off between multiple reduction and 
reflector distortion must be determined by trial-and-error 
for particular applications. 

4. Multiple Distortion 

Since homomorphic dereverberation is theoretically based 
on a model of strictly impulsive multiples, it is important to 
observe the performance of this technique in the more realistic 



Figure 49 





(b) 




(c) 
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0.5 second section of cepstra at multiple location, 
(a) MSR=20 (b) MSR=1 (c) MSR-.06 
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case of an extended bottom interaction mechanism. We employ 
the same bottom response functional form discussed earlier 
in this chapter. The seismogram of figure 50a contains 
distorted multiples due to a bottom response of .045 seconds 
duration. Reflectors are present at 2.7 and 3.5 seconds. 

Figure 50b shows the region of the cepstrxim centered at the 
first multiple location. The peaks are much broader than those 
of figure 49. The dominant energy is concentrated near the 
multiple location as before but it now appears to be more 
distributed. 

Table 1 summarizes the effects of varying the stopband 
width and location in the cepstrum of figure 50b. First and 
third multiple energy removed varies widely while the second 
multiple is not significantly affected. by the passband dimensions. 
In some cases, extending the stopband decreases the amount of 
multiple energy removed. Zeroing the .94-. 98 second region 
results in greatly increased reduction of the third multiple 
(74-78%) but first multiple reduction is degraded by 10 - 12 %. 

This behavior, which has been observed in several cases, is 
not predicted by the theory and is thought to be a computational 
effect. The high amplitude oscillations which are dominant in 
the left side (.75-1.0 second) of figure 50b, or the low 
frequency "drift" which is apparent throughout the section may 
)30 computational noise which contributes to this phenomenon. 

Several observations can be made in this case. Although 
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Figure 50 (a) Seismogram with distorted multiples 

(b) 0.5 second section of cepstrum centered at 

first multiple location. 
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stopband Limits ENERGY REMOVED 

( S0 c ) 

1st MULT 2nd MULT 3rd MULT 1st REEL. 2nd REEL. 















. 94-1.06 


.75 


.48 


.78 


-.13 


-.01 


.96-1.06 


.78 


.44 


.74 


-.13 


.01 


.96-1.00 


.65 


.49 


.56 


-.12 


-.03 


.98-1.02 


.79 


.44 


.45 


-.12 


-.03 


.98-1.04 


.88 


.44 


.43 


-.12 


-.03 


.98-1.06 


.89 


.45 ! 


.54 


-.12 


-.02 


.99-1.04 


.78 


.42 


.39 


-.12 


-.05 


1.0-1.04 


.64 


.45 


.48 


-.13 


-.06 


1.02-1.06 


.50 


.45 


.50 


1 -.13 

1 


-.06 ! 

1 















TABLE 1 

EEEECTS OE STOPBAND LIMITS ON PEREORMAfJCE 
EOR THE SIGNAL OE EIGUP£ 50. 
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performance is somewhat erratic^ a stopband of 40-80 msec 
which straddles the first multiple onset time yields significant 
multiple reduction. The peak performance in this case is about 
10% below that achieved v/ith undistorted multiples but the 
sensitivity to minor stopband variations is less dramatic. 
Reflector distortion does not increase significantly due to the 
presence of distorted multiples. 

Figure 51 shov;s processed seismograms resulting from the 
application of various stopbands to the cepstrum of figure 50. 

Finally, we note that the effects of aperiodic! ty have 
not been discussed in the distorted multiple case. Due to 
limitations of the earth model it was not possible to investigate 
these effects. Thusfar, however, we have relaxed the theoretical 
assumptions regarding periodicity and coherence individually, 
and found that homomorphic de reverberation is not critically 
sensitive in either case. We surmise that the combination of 
these effects would not be catastrophic to performance. Further 
investigation is merited since this topic has an important 
bearing on the effectiveness of the homomorphic technique in 
shallow water dereverberation where later multiples are 
significant. 

5. Water Travel Time Estimate 

In practice, passband location must be determined by 
estimation of water column travel tim.e. It is apparent from 
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Figure 51 Processed results of the seismogram in figure 50a 

for three different cepstral stopbands. (a) 1.02-1.06 sec 
(b) 1.0-1.04 sec (c) .98-1.06 sec 
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the preceding discussion that the importance of accurate 
bottom estimation depends on the seismic environment. In 
particular, we have seen that MSR and bottom characteristics seen 
to affect cepstral energy distribution. Proximity of multiples 
to important reflectors also dictates resolution constraints 
on travel time estimation. For the cases analyzed here the 
required travel time estim.ation accuracy would be 20-40 msec 
since, as determined in the foregoing discussion, a stopband 
of 40-80 msec is generally required for effective dereverberation 
Travel time estimation error of more than half the stopband 
width causes the major multiple contributions to be excluded 
from the stopband. This degree of accuracy can easily be 
attained with existing bottom tracking algorithms. 

The reflector/multiple resolution -implied by the required 
stopband widths is also about 40-80 msec for the data tested. 

A discussion and several examples of reflector/multiple resolu- 
tion are included in the follov;ing section. 

6 Reflector/Multiple Overlap 

The follov/ing six figures illustrate homomorphic dereverbera 
tion of signals in which reflectors and multiples .are closely 
situated . 

Figure 52a shows a seismogram with reflectors at 1.6, 2.02 
and 3.0 seconds. The first and second multiples directly inter- 
fere with reflectors and the third multiple is barely visible 
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Figure 52 (a) Unprocessed seismogram 

(b) Result of cepstral notch filtering; 
stopband 

(c) .98-1.06 sec stopband 



.96-1.02 sec 
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at 4.0 seconds. The result of applying a cepstral stopband at 
.96-1.02 seconds is seen in figure 52b. Most of the energy 
near two seconds has been removed but a small reflection is 
still visible. The signal of figure 52c results from a stop- 
band of .98-1.06 seconds which spans .the time domain reflector 
location as well as the multiple location. Residual energy 
is still present after this processing which indicates that 
some of the reflector and first multiple cepstrum contributions 
are distributed beyond the immediate vicinity of the time 
domain locations . 

The seismogram of figure 53a contains reflectors at 1.6, 
1.95 and 3.0 seconds. Each reflector is clearly evident after 
processing C. 9 8-1. 06 second stopband) and first multiple energy 
has been greatly reduced as shown in figure 53b. Some of the 
removed energy may be due to the 1.95 second reflector; however, 
the first multiple in this example is very coherent and its 
energy is more likely to be localized near 1.0 second in the 
cepstrum. Extension of the stopband closer to the reflector 
location C. 96-1. 06 second) causes complete removal of the 
reflector as seen in figure 52c. Visible reduction of the 
second multiple at 3.0 seconds and an internal multiple at 2.7 
seconds is apparent in figure 52b. 

The seismogram of figure 54a has the same reflector place- 
ment as that of figure 53a but the multiples in this case are 
smaller and less coherent. Application of a cepstral stopband 
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Figure 53 (a) Unprocessed seisnogram (b) -Result of cepstral 

notch filtering; stopband .98-1.06 sec. (c) .96-1.06 

sec stopband. 
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3.0 4.0 5.0 



(a) Unproccscecl seismogram 

(b) Result of cepstral notch filtering; stopband 
1.0-1.06 sec. 



Figure 54 
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at 1.00-1.06 second yields (figure 54b) visible reduction 
of all three multiples and well resolved reflections at the 
proper locations . 

Figure 55a illustrates a seismogram with considerable 
multiple distortion and reflectors at 1.5, 2.1 and 3.0 seconds. 

A wide cepstral stopband, .98-1.08 second, has a very minor 
effect on the first multiple as seen in figure 55b. Extension 
of the stopband to 1.12 seconds causes removal of virtually 
all signal energy in the region. The relatively wide (.1 second) 
ref lector/multiple separation cannot be effectively exploited 
in this case because the very incoherent multiple in a low MSR 
signal requires a large stopband for effective removal (see 
figure 48) . 

The result of longpass filtering the cepstrum of the signal 
in figure 55a is shown in figure 56. All cepstral contributions 
prior to 1.02 seconds have been set to zero. The bottom and 
later reflectors are clearly visible but the 1.5 second reflector 
has been rem.oved. This illustrates one drav;back to longpass 
dereverberation in deeper v;ater. This effect may be acceptable 
in some situations, however. For instance, good resolution of 
closely spaced, deep reflectors may be obtained by longpass 
filtering whereas the earlier reflectors are frequently obvious 
before processing. 

Figure 57 shows longpass results for a signal with impulsive 
multiples and very sharp reflectors at 2.2 and 3.0 seconds. 
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Figure 55 (a) Unprocef=;sed r.eisrnogram. 

Cb) Result of cepstral notch filtering; stopband 
.98-1.08 sec 
(c) .98-1.12 sec 
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Figure 56 (a) Unprocessed seismogram 

(b) Result of longpass filtering cepstrum; cutoff 
quefrency 1.02 sec. 
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Unprocessed seisnograin 

Result of longpass filtering; cutoff quefrency 
1.01 sec. 



Figure 57 (a) 

(b) 
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The stopband limit is 1.01 seconds in this case. Multiple 
removal is complete and both reflectors are clearly visible. 
Longpass filtering was found to be most effective for signals 
of this type. Very sharp reflectors are more clearly visible 
than more distorted reflectors of comparable energy, especially 
in noisy signals. Longpass filtering of noisy signals is 
discussed later in this chapter. 

The low frequency noise present in the reflector region 
of figure 57 is typical of that observed in several longpass 
results. The cepstra of figures 49 and 50 contain similar 
components. No attempt was made in this analysis to remove 
this type of noise. The reason for its presence has not been 
determined . 

7. Additive White Noise and Filtering 

The effects of additive noise are not addressed in the 
formulation of the homomorphic system for deconvolution since 
there are no apparent cepstral properties of noise which can be 
exploited for its removal. Linear filtering is a more suitable 
way of reducing additive noise in individual signals. This can 
be performed prior to, or in conjunction with, homomorphic 
deconvolution. The effects of this combined processing in the 
special case of Gaussian white noise have been investigated and 
are discussed here. The data presented represent a relatively 
small number of experiments performed with synthetic data and 
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noise generated as previously discussed. These results cannot 
be generalized categorically due to the relatively narrow 
scope of the experiments and the lack, of precise theoretical 
characterization of noise properties under homomorphic trans- 
formation. The results do indicate perfoimiance trends and 
computational effects due to additive noise and filtering. 

Since noisy signals usually undergo linear filtering prior 
to dereverberation, we begin by discussing an important 
computational issue which arises when signals are bandpass 
filtered. Such filtering introduces spectral regions (stopbands) 
containing little or no spectral energy. Recall that the homo- 
morphic transformation involves computing the logarithm of the 
Fourier transform of the signal. Since the logarithm of zero 
is not defined, this computation is not possible in frequency 
bands where the Fourier transform is zero. In the case of 
lowpass filtered signals this problem can frequently be overcome 
by resampling after filtering. Figures 58 a, b and c illustrate 
this process. If the sampling frequency is decreased to the 
Nyquist rate implied by the filter cutoff frequency, the 
resulting discrete Fourier transform will- include only the 
frequency components in the passband region. Figure 58d shows 
a distribution of spectral energy which is not readily amenable 
to elimination of the zero region. In this case tha baseband 
j^ggion is zero so that resampling would not be effective. 
Investigation of this problem is currently in progress 
(Tribolet [11] ) . 
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Figure 58 (a) Unfiltered spectrum sampled at the Nyquist rate. 

(b) Spectrum of (a) after lowpass filtering at 
without resampling. 

(c) Filtered spectrum after reampling. 

(d) Spectrum with zeros in the baseband region. 
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It was found in this analysis that computing the cepstrum 
of an over-sampled signal is an unstable procedure which, at best, 
requires extensive computation time. The low energy regions 
of the spectrum lead to spurious phase derivative values . 

Hence, the phase unwrapping algorithm proceeds at small step 
sizes, requiring computation of many intermediate values of the 
discrete Fourier transform. In many cases the computer algorithm 
could not produce an unwrapped phase which was in acceptable 
agreement w’ith the principal value. Cepstra of over-sampled 
signals which were successfully computed, however, generally 
yielded dereverberation results comparable to those of properly 
sampled signals (see table 2) . 

Addition of white noise was found to have a computational 
effect similar to that described above; phase unwrapping time 
is significantly increased. Heavy weighting (w =.98-. 99) 
was found to reduce computation time considerably for noisy 
signals, although some signals require small phase integration 
step sizes in isolated sections even v;hen substantial weighting 
is applied. Recall from Chapter II that there is no assurance 
that the log-spectrum is adequately sampled, even after lowpass 
filtering of the signal. Hence, the integration of the phase 
derivative cannot be expected to proceed quickly in all cases. 

Smoothing of the phase derivative was attempted to 
compensate for the effects of noise. A three-point moving 
average was applied prior to integration. The resulting inverse 
cepstra bore no resemblance to the original seismograms. 
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SAMPLING 

INTERVAL 

(MSEC) 



FILTER 
CUTOFF (Hz) 



FRACTION OF ENERGY REMOVED 
1st MULT 2nd MULT 3rd MULT 













4.88 


50 


.91 


.064 


.25 


9.77 


50 


-.02 


.12 


.24 


19.53 


50 


.94 


.044 


.25 


4.88 


20 


.94 


.083 


.15 


9.77 


20 


. 93 


.48 

( 


.13 


19.53 


20 


.98 ! .12 


.55 













TABLE 2 



EFFECT OF RESAMPLING ON MULTIPLE REMOVAL FOR A 
NOISELESS SEISMOGRAI4 LOWPASS FILTERED AT 50 and 
20 Hz. DEREVERBERATION WAS ACCOMPLISHED BY APPLY- 
ING AN 80 MSEC CEPSTRAL STOPBAI^D AT THE FIRST 
MULTIPLE LOCATION. 
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Figure 59 summarizes the observed effects of noise on 
fij^st multiple removal by homomorphic processing. Each 
seismogram was lowpass filtered at 50 Hz and resampled by a 
factor of four prior to multiple removal. Percentage of 
niultiple energy removed shows a consistent decrease with increas- 
ing noise level . The rates of decrease and amounts of energy 
removed are seen to vary widely from one signal to another. 
Examples of noisy signals before and after processing are shown 
in figures 60 and 61. Significant reduction of the first multiple 
is evident in both examples. In the first case, the noise level 
is moderate and multiple reduction has a marked effect on visual 
quality of the signal. The noise level in figure 61 is consider- 
ably higher, resulting in marginal improvement due to dereverbera- 
tion. 

Second and third multiple energy removed was found to 
decrease generally with decreasing noise also, and in some 
high noise cases the later multiples were actually enhanced. 

Data for a typical signal are tabulated below. 



STANDARD 
DEVIATION 
OF NOISE 



ENERGY REMOVED 
2nd MULT. 3rd MULT. 









0 


.29 


.32 


25 


.19 


.15 


50 


.04 


.02 


100 


1—1 

<N 

• 

1 


-.17 









TABLE 3 
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STANDARD DEVIATION OF ADDED NOISE 



Figure 59 Effect of added noise on homoiuorphic dereverberation. 
All signals are lowpass filtered at 50 Hz ^nd 
resampled. 
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1.0 2.0 3.0 4,0 5.0 



Figure 60 (a) Seismogram with moderately high noise, F =50 Hz. 

(b) Result of cepstral notch filtering (.8-l54 sec). 





Figure 61 Seisnogram with very high noise, F =50 Hz. 

(b) Result of cepstral notch filtering^ (. 80-1 . 04 sec). 
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Reference to figure 61 reveals that for the higher noise 
amplitudes the later multiples are buried so that their removal 
is not important. 

Reflector distortion/ summarized for a typical case in 
table 4 , was found to be generally more severe than in noise- 
less signals but rarely more than 25%. 



STANDARD 

DEVIATION REFLECTOR ENERGY REMOVED 

OF NOISE 











0 


-.18 


.02 


O 

• 


25 


-.22 


• -.25 


-.22 


50 


1 

• 

o 

L 


-.14 


-.05 


100 


.14 


-.12 


o 

• 

1 











TABLE 4 

EFFECT OF NOISE ON REFLECTOR DISTORTION 
FOR A TYPICAL SEISMOGRAM WITH 3 REFLECTORS, 
FILTERED AT 50 Hz. 



The above results indicate that the addition of noise 
leads to decreasing performance with respect to all three 
criteria. The behavior is somewhat spurious and frequently 
exhibits wide variation from signal to signal. 
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Several lowpass filter bandwidths (20, 30, 50, 70, and 
90 Hz) were applied to the signals evaluated in figure 59. 

The results are shown in figure 62, In each case, first multiple 
removal is least effective when the signals are prefiltered at 
50 Hz. Figures 63 and 64 illustrate this behavior. The same 
seismogram is shown in figures 63a and 64a, lowpass filtered 
at 50 and 20 Hz respectively. The difference in appearance is 
dramatic. Resolution of the reflectors on either side of the 
first multiple (at 1.55 and 2.4 sec) is greatly improved by 
dereverberation in the latter case while figure 63 shows very 
little visual improvement. 

The behavior illustrated in these figures cannot be fully 
explained on the basis of the data available. The signals 
tested have dissimilar reflector locations and varying amounts 
of multiple distortion which implies that the similar performance 
dips are not due to similarities in signal configuration. The 
observed behavior may be due to the properties of the source 
signature (which is common to all three seismograms) or the 
characteristics of the recursive third order Butterworth 
filter employed. More data using different filter routines 
and a wide range of signal characteristics is needed to 
completely characterize this behavior. The results obtained 
here suggest that a filter bandwidth very close to the bandwidth 
of the noiseless signal leads to the best dereverberation 
performance. The same result was obtained for the TDL algorithm. 



PERCENT OF FIRST MULTIPLE ENERGY REMOVED 
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FILTER CUTOFF FREQUENCY (Hz) 



Figure 62 Effects of lov;pass filtering signals with added 

white noise prior to homomorphic dereverberation . 
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(a) Noisy seismogram lowpass filtered at 50 Hz. and 
resampled at 51 Kz. 

Cb) Result of cepstral notch filtering; stopband 
.8-1.04 sec. 



Figure 63 
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Figure 64 Noisy seismogram lowpass filtered at 20 Hz and 

resampled at 51 Hz. Cb) Result of cepstral filtering; 
stopband .8-1.04 sec. 
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Additive noise was found to have a particularly adverse 
effect on longpass filtered signals. Noise, which may be 
relatively low in the received signal, tends to be amplified 
with respect to reflectors in the process of longpass filtering, 
especially in the later portion of the signal. Although long- 
pass filtering removes a large part of the noise energy with 
the source, the overall effect is usually a decrease in SNR. 

The advantage of longpass filtering, which includes source 
deconvolution, is that greater resolution of close reflectors 
can be achieved. It was found that this approach is worthwhile 
in low noise signals but not effective when the white noise 
level is significant with respect to reflector amplitudes. 

Proper resampling after prefiltering was found to be very 
important for successful longpass filtering. Figure 65a 
illustrates the effect of longpass filtering the cepstrum of a 
noisy signal which was first lowpass filtered at 20 Hz but 
not resampled. The filtered signal contains relatively low 
noise and the sampling frequency is 205 Hz. The processed 
result of figure 65b is useless due to the high sampling rate. 
Reduction of the sampling frequency to 102.5 Hz leads to the 
processed signal of figure 66a. The 3.5 second reflector is 
clear but the high background noise almost obscures the 2.7 
second reflector. Multiple removal is complete. Resampling 
to 51 Hz, which is approximately the Nyquist rate in this case, 
leads to some improvement (figure 66b) but the SNR is much lower 
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Figure 65 (a) Noisy seismogram lowpass filtered at 20 Hz but 

not resampled. 

(b) Result of longpass filtering the cepstrum of (a) . 
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Figure 66 Longpass processing results Ca) for the signal of 

figure 65a, resampled at 102.5 Hz before processing, 
(b) for the signal of figure 65a, resampled at 51 Kz 
before processing (c) for a signal identical to 
figure 65a with one half the v;hite noise level, 
resampled at 51 Hz. 
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than before dereverberation (figure 65a) , Figure 66c results 
from processing identical to that of figure 66b, on the same 
signal, with the noise amplitude halved. Both reflectors are 
clear and the multiples have been removed but the background 
noise is much higher even than in the signal of figure 65a. 

Thus, longpass filtering of noisy signals is seen to involve 
a trade-off between effective dereverberation and decrease in 
SNR. For signals with moderate to heavy noise the reduction in 
SNR was found to be unacceptable in the cases tested. 

Although this subject was not extensively investigated 
there is some indication that lowpass filtering can be employed 
to improve the results of subsequent longpass processing. 

Figures 67 and 68 illustrate the longpass processing of a 
noisy signal after lowpass prefiltering at (a) 70 Hz, (b) 50 Hz 
and (c) 30 Hz. In each of the signals in figure 68 the multiples 
have been removed but the reflector resolution improves 
considerably from (a) to (c) . These results are reasonable 
in that improvement of performance coincides with increasing 
rejection of out-of-band noise; however, the filter and signal 
characteristics must be studied more closely to explain this 
behavior accurately. 

D. Comparative Examples of Processing Results 

The foregoing results and discussion illustrate the 
performance of the TDL and homomorphic dereverberation algorithms 
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Figure 67 Noisy sei sinogram lowpass filtered at three different 
frequencies (a) 70 Kz (b) 50 Kz (c) 30 Kz. Each has 
been resampled at 51 Hz which is the approximate 
Nyquist rate for the signal without noise. 
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Figure 68 Results of longpass filtering the seismograms of 
figure 65 a, b and c, respectively. 
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for a variety of processing conditions. Several trends in 
performance are apparent. Before proceeding to a summary 
and discussion of the relative strengths and weaknesses, we 
present several examples which allow direct visual comparison 
of the techniques. Each of the following figures contains 
(a) an unprocessed seismogram and the two processed results 
obtained by (b) homomorphic and (c) TDL filtering. 

Figure 69 illustrates a noiseless seismogram with multiple/ 
reflector overlap at both 2 and 3 seconds. Both algorithms 
leave a small amount of energy at the first location and effect 
only slight reduction at the second. In general, both methods 
were found to eradicate reflectors which are extremely close 
to the first multiple and retain signal components which closely 
coincide with later multiples. 

In figure 70 the multiple onset occurs .1 second before 
the reflector and the MSR is considerably lower than in the 
previous figure. In this case the separation is great enough 
that both methods retain a significant amount of signal energy 
near 2.1 seconds. The homomorphic result is considerably 
sharper although very little energy has been removed. Multiple/ 
reflector separation of .1 second was found to be the approxi- 
mate resolution limit of both techniques when reflector onset 
is later than multiple onset and travel time is estimated 
accurately . 

Figure 71 shows a reflector at 1.95 seconds, slightly 
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(a) Unprocessed seismogram with reflectors at 1.5, 2.02 
and 3.0 seconds. 

(b) Result of homomorphic processing. 

(c) Result of TDL processing. 



Figure 69 
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Figure 70 (a) Unprocessed seismogram with reflectors at 1.5, 2.1 

and 3.0 sec. 

(b) Result of homomorphic processing 
Cc) Result of TDL processing. 
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Figure 71 



Ca) Unprocessed seismogram with reflectors at 1.6, 1.95 
and 3.0 sec. 

(b) Result of homomorphic processing. 

(c) Result of TDL processing. 
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before the multiple. Effective dereverberation is accomplished 
in both b and c. As in the previous figure the reflector is 
better resolved by homomorphic processing. Resolution of both 
techniques was generally observed to be slightly better when 
reflector onset is earlier than multiple onset, provided travel 
time is estimated accurately. In such cases of severe inter- 
ference the homomorphic filtering generally yields more distinct 
reflections. A further example of this behavior is shown in 
figure 72. 

The following three figures illustrate dereverberation of 
noisy signals. Each seismogram has been lowpass filtered at 
30 Hz and the homomorphic outputs as shown have been resampled 
at 51 Hz. In this first case (figure 73) the performance of 
both methods is comparable. The first multiple is almost 
completely removed while other regions of the signal are not 
visibly affected. Figure 74 also shows comparable multiple 
removal, however, the resolution of the reflectors at 2.6 and 
3.4 seconds is somewhat better after TDL filtering. The 
homomorphic algorithm was found to produce higher random noise 
spikes than the TDL filter. This effect is present in figure 
74 and again in figure 75. In both examples the homomorphic 
method achieves slightly better multiple removal but the 
overall noise level in the result appears higher. 
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(a) Unprocessed ooismogram v;ith reflectors at 1.6 
1.95 and 3.0 sec. 

Cb) Result of homomorphic processing. 

Cc) Result of TDL processing. 



Figure 72 
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(a) 






1.0 2.0 3.0 4.0 5.0 



Figure 73 Ca) Noisy seismogram with reflectors at 2.7 and 

3.5 sec, and lo^vpass filtered at 30 Hz. 

(b) Result of homomorphic processing. 

(c) Result of TDL processing. 
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Figure 74 Ca) Noisy seismogram v/ith reflectors at 1.55, 2.4 and 

3.25 sec, lowpass filtered at 30 Hz. 

Cb) Result of homomorphic processing. 

(c) Result of TDL processing. 
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(c) 



1.0 2.0 3.0 4.0 5.0 

Ca) Noisy seisrnograra vrith reflectors at 2.6 and 3.4 
sec, lov/pass filtered at 30 Hz. 

(b) Result of homomorphic processing. 

(c) Result of TDL processing. 



Figure 75 
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CHAPTER V 

DISCUSSION AND SUMMARY 



In the preceding chapters we have (1) examined the 
theoretical structures of the TDL and homomorphic dereverbera- 
tion techniques, (2) established a comprehensive set of perfor- 
mance criteria, and (3) presented the results of applying both 
methods to synthetic data. Our approach has been essentially 
that of perturbation analysis. Through variation of environ- 
mental and signal processing parameters we have observed 
performance trends due to deviations from the ideal theoretical 
models upon which the methods are based. In a more practical 
sense the parameter variations simulate a range of seismic 
processing conditions. Since there are many different environ- 
ments in which these algorithms may be applied, we have not 
emphasized the absolute performance figures obtained from the 
simple, synthetic data utilized here. Rather, we have tried 
to present a behavior profile of both algorithms which indicates 
the basic trends and sensitivities with respect to a number 
of parameters, interpreted in terms of their theoretical 
structures and assumptions. This approach is intended to give 
a more general indication of the dereverberation performance 
to be expected in different situations. We conclude with a 
summary and discussion of the comparative results presented in 
Chapter IV. 
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In low noise, deep water seismograms the methods have 
been found to be comparable for reducing dominant first 
multiples (by 75-95% in most cases) . The TDL filter requires 
accurate signal statistics including an estimate of the multiple- 
reflector crosscorrelation function, which must be approximated. 
This leads to degraded performance in shallow water situations 
where significant reflector energy is within the crosscorrelation 
window. The homomorphic method requires no statistical character- 
ization of the signal and thus has no similar performance 
degradation in shallow water; however, the three-stage cepstral 
transformation requires extensive computation which may be an 
important limitation for at-sea processing systems. (This 
issue v/ill be discussed in more 'detail later.) 

Although the effects of aperiodicity could not be thoroughly 
evaluated experimentally the derived result expressed in 
equation (11) suggests that the homomorphic algorithm has the 
potential to reduce later, aperiodic multiples. The combina- 
tion of such processing with source deconvolution by longpass 
filtering the cepstra of shallow water seismograms appears to 
be the most promising application of homomorphic dereverbera- 
tion. More extensive practical evaluation is needed in this 
area. 

Closely spaced, aperiodic multiples destroy the coherence 
of the approximated crosscorrelation function at shifts near 
the two-way travel time which, again, limits the effectiveness 
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of the TDL filter in shallow water signals. 

Increasing white noise level causes monotonically decreas- 
ing performance in both methods as illsustrated in figures 38 
and 57 . The TDL performance fell off more slowly with noise 
when both techniques were applied to. similar signals. Bandpass 
filtering leads to a consistently higher percentage of multiple 
reduction by TDL processing of noisy signals. Filter effects 
on homomorphic processing are more complex. The cases evaluated 
indicate that multiple energy removed is not a monotonic func- 
tion of filter cutoff frequency. Considerably more data are 
required to determine the precise effects of filter bandwidth 
and phase characteristics. The resampling which was found to 
be helpful after bandpass filtering may have a .detrimental 
effect on visual record quality, so that interpolation may be 
desirable in some cases. 

Reflector distortion does not appear to be a problem for 
either technique except in cases where overlap is severe. In 
most cases tested less than 10% of the reflector energy was 
removed. Close proximity of reflectors to the first multiple 
frequently leads to severe distortion by. both methods due to 
the resulting bias of the crosscorrelation function and the 
lack of sufficient cepstral separation. Reflectors close to 
later multiples are usually well preserved. This behavior 
suggests two applications of the homomorphic algorithm. First, 
when regions of geological interest occur after the onset of 
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the first multiple, a wide pepstral stopband can be employed at 
the first multiple location to reduce later multiples without 
distorting reflectors. In very shallow water the stopband may 
be extended to include more than one multiple. A second 
possibility for avoiding reflector distortion is the use of 
weighting coefficients greater than 1.0 to exploit the properties 
of mixed phase sequences. The object of this weighting is to 
make the reflector train mixed phase while keeping the multiple 
sequence minimum phase. This appears to be feasible in many 
situations since the z-plane zero of the multiple sequence is 
usually well inside the unit circle. Moving some of the reflector 
train zeros outside of the unit circle (i.e., making it mixed 
phase) will, in general, cause some of the cepstral energy due 
to the reflectors to occupy the negative quefrency region. 

Even if the reflector train has a maximum phase component 
before weighting the same effect can be expected. Thus, the 
amount of reflector energy near the first multiple location 
may be reduced. Although there is no guarantee that the 
resulting notch filtered cepstrum will transform to a seismogram 
with less distortion, this technique appears to be worthy of 
investigation . 

We recall one other reflector distortion effect which was 
seen in Chapter IV. We saw in figure 55 that dereverberation 
by longpass filtering completely removes reflectors occurring 
prior to multiple onset. It was noted that this effect may be 
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acceptable if it leads to better resolution of smaller, deep 
reflectors . 

In terms of our third criterion, visual improvement of the 
signal, both techniques were seen to have advantages and 
disadvantages. Homomorphic processing usually resulted in 
better resolution of interfering signal components in the lov; 
noise signals processed. Increasing noise, however, was seen 
to cause homomorphic results more prone to random noise spikes 
which degrade the interpretability of the record. The results 
of TDL filtering had somewhat better visual resolution in the 
noisier signals processed. As noted previously the visual 
advantage of the TDL in this case is partially due to the 
noisier appearance of the resampled signals produced by 
homomorphic processing. 

Longpass filtering was seen to provide the most effective 
dereverberation, the best reflector resolution and the best 
overall visual quality in ideal cases. Unfortunately it 
degenerates quickly with noise and could not be successfully 
applied to very noisy signals or signals with important geological 
regions preceding the multiple onset. Further research and 
experience may well lead to more extensive applicability of 
this technique. 

The relative simplicity of the TDL algorithm makes it 
much more desirable from a computational standpoint. TDL 
dereverberation of a 1024 point signal can be accomplished in 
3 seconds or less for the operator lengths typically required. 
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The three major computational steps are the correlation 
operations, solution of the estimator equations and convolu- 
tion of the operator with the signal. 

The homomorphic computations include weighting, four FFT's 
computation of the complex logarithm, phase unwrapping, linear 
filtering, complex exponentiation and unweighting. This 
algorithm can be expected to take 20 seconds or more for a 
1024 point sequence on a small processing computer. In this 
analysis the phase unwrapping computation took over one minute 
for some noisy signals. These figures are highly dependent 
on hardware available and programming efficiency but, in 
general, homom.orphic dereverberation is several times slower 
than the TDL algorithm. Special purpose hardv/are could be 
used to reduce homomorphic computation time significantly, 
but the method has not been implemented for processing on a 
large scale thusfar. 

Storage requirements for the homomorphic algorithm vary 
with the FFT routine used, method of cepstrum computation and 
cepstrum length. The program used for this analysis requires 
about 12 * N bytes of core and 4 * N bytes of disc storage, 
where N is the cepstrum length. N was twice the signal length 
in the cepstra computed. Shorter cepstrum lengths, as 
determined by trial-and-error , may produce results with 
acceptably lov; aliasing in many cases. The TDL dereverberation 
program used requires about 6 * L bytes of core, where L is 
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the data sequence length. No disc storage is required in 
this computation. These storage requirements apply to a float- 
ing point processing scheme on a machine (HP-2100) which uses 
four byte floating point words. 

In conclusion, we make some general observations concern- 
ing the results of this analysis. 

The TDL dereverberation scheme is a simple and efficient 
technique which has demonstrated effectiveness in removing 
deep water multiples. The analytical structure is well under- 
stood and its performance characteristics have been explained 
here in terms of that structure. Further refinements in its 
implementation may be possible but its potential is essentially 
clear at this point. 

Homomorphic dereverberation is complex, relatively untested 
and requires extensive computation. It has been shown here 
to be effective on synthetic data. It appears to be particularly 
promising for shallow water dereverberation. The complexity 
of the method leads to a number of possible analytic and 
computational techniques which can be utilized in its applica- 
tion. Certainly, its full potential has not yet been determined 
and further investigation is warranted. 



REFERENCES 



(1) Backus, M. , "Water Reverberations - Their Nature and 
Elimination," Geophysics, v.24, no. 2 . pp. 233-261 (April. 
1959) . 

(2) Mayne, W., "Common Reflection Point Horizontal Data 
Stacking Technique," Geophysics, v.27, no. 6, pp. 927-938 
(December, 1962) . 

(3) Baggeroer, A. B., "Tapped Delay Line Models for Dereverbera- 
tion of Deep Water Multiples," IEEE Trans, on Geophys . Elec. 
V. G.E. 12, no. 2, pp. 33-54 (April 1974). 

(4) Van Trees, H. L. , Detection, Estimation and Modulation 
Theory, Part III , John Wiley & Sons, Inc.: New York, 1971. 

(5) Stoffa, P. W. , Buhl, P., Bryan, E.M., "The Application of 
Homomorphic Deconvolution to Shallow Water Marine Seismology 
Geophysics , v.39, no. 4, pp. 401-426 (August, 1974). 

(6) Oppenheim, A. V., Schafer, R. W. , Digital Signal Processing , 
Prentice Hall, Inc.: Englewood Cliffs, N.J., 1975. 

(7) Schafer, R. W. , "Echo Removal by Discrete Generalized Linear 
Filtering," Research Lab. of Electronics, M.I.T., Tech. Rep. 
466 (February, 1969) . 

(8) Ulrych, T. J. , "Application of Homomorphic Deconvolution 
to Seismology," Geophysics , v.36, no. 4, pp. 650-660 
(August, 1971) . 

(9) Tribolet, J. M. , "A Nev; Phase Unv;rapping Algorithm," 
submitted to IEEE Trans, on Acoustics, Speech and Signal 
Processing , (March, 1976). 

(10) Theriault, K. B., "Optimum Arrival-Time Estimation in 
Exploration Seismology," Ocean '74 IEEE International 
Conference on Engineering in the Ocean Environment , IEEE , 
Inc. , New York (1974) . 



(11) Tribolet, J. M. , personal communication, April 1976. 



-163- 



APPENDIX A 

COMPUTATION OF THE PHASE DERIVATIVE OF THE Z-TRANSFORI4 



When computing the complex cepstrum of a sequence, x(n), 
it is necessary to determine the unique, continuous phase of 
X(z) . One v/ay of obtaining the continuous phase is to first 
compute its derivative and then integrate numerigally. The 
computation of the phase derivative from x(n) is discussed 
in detail here. 

We begin v;ith the z-transform of x(n), 

00 

X(z) = I x(n) z ^ = Xj^(z) + jXj(z). 

n=-oo 

Taking the complex natural logarithm 

/s 

log X(z) H x(z) = log lx(z)| + j arg X(z) 



We see that the phase of X(z) is equal to the imaginary 
part of its natural logarithm. The derivative of X^(z) can be 
expressed in terms of easily computable quantities. 



dX (z) 
dz 



log[X(z)] 

dz 



X' (z ) 
X(z) 



( 1 ) 



where the prime indicates differentiation with respect to z. 

^ I 

Expanding (1) and solving for X^(z), 
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X^(z) + jXj(z) = 



X^(z) + jxj(z) 

Xj^(z) + jXj(z) 



X*(z) - jX (z) . 

X^(z) = -I + j X^(z) 

Xj^(z) + jXj(z) 



Separating the lyiS into real and iinaginary parts, 

X (z) X'(z) - X^(z) x'(z) 

Xj(z) = ± ± 5_ 

Xj^^(z) + Xj(z) 



+ j 



Xr(z) 



(Xr(z) Xr(z) + Xj(z) Xj(z)) 
Xr(z) + Xj(z) 



The real part yields an expression for the phase derivative, 



X (z) X|(z) - X^(z) X„(z) 

xj(z) = i ± 

X^(z) + X^(z) 



( 2 ) 



Since the z-transforra is actually evaluated on the unit circle 
using the discrete Fourier transform (DFT) we set z = e^'^. 

X fe^“) x;(e3‘^) - X,(e^‘^) x'(e^“) 



xj(e^“) = — 



R 



X^(e^“) + 



( 3 ) 



Derivatives with respect to e^'^ may be replaced by ^ 



Since 
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dX(e^“) 


j(0 

— 1 /-N -f 


^dX(ei“)' 


d(o 

1 J 


- J e 


de^“ 

J 



and we thus have a common factor of j e^^ in each term of (3) . 

Hence, we need only compute the real and imaginary parts 
of X(e^^) and ^ (x(e^^)} and combine them as indicated in (3) 
The derivative of X(e^^) is easily obtained from the sequence 
n x(n) as follov;s: 

00 n to 

V / jw» V , > -jojn X(e-^ ) 

X (e-* ) = ^ n x(n) e -• = ]d 

n=-°o 



Re[X^(e^“)] = 


-x;(e^“) 


Im[X^(e^“)] = 


x^(e^“). 



The required computation in terms of the DFT is then 



xj (k) 



X^(k) + Xj (k) 
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